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Streszczenie 



The aim of the thesis is to present and analyze two particular problems of 
transport in porous media flow. The first of them is related to the process of sa- 
turation of porous building materials. Recently, M. Kiintz and E Lavalee, using a 
computer model of this process, have concluded that the anomalous diffusion 
assumption is correct. In this thesis I present an alternative explanation of this 
results without any refer to anomalous diffusion. The second part of the thesis 
covers the numerical analysis of the tortuosity of the flow - one of a very intere- 
sting physical macroscopic variables characterizing transport in porous media. 

Celem niniejszej rozprawy jest przedstawienie oraz analiza dwoch szczego- 
lowych problemow zwiazanych z transportem plynow w osrodku porowatym. 
Pierwszy z nich zwiazany jest z zagadnieniem zwilzania materialow budowla- 
nych. Niedawno M. Kiintz i P. Lavalee opracowali model komputerowy tego pro- 
cesu i stwierdzili, ze mamy w nim do czynienia z dyfuzja. anomaln^. W niniejszej 
rozprawie przedstawiam alternatywne wyjasnienie tych wynikow bez odwofy- 
wania sie do dyfuzji anomalnej. Drugie zagadnienie dotyczy szczegolowej ana- 
lizy numerycznej kretosci przeplywu - jednej z ciekawszych makroskopowych 
wielkosci fizycznych charakteryzujacych transport wosrodkach porowatych. 
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Rozdziaf 
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Wst§p 



1.1. Wprowadzenie 

Transport ciepla, cieczy i gazow przez osrodki porowate jest zjawiskiem wszech- 
obecnym, a badanie jego wtasciwosci ma ogromne znaczenie dla wielu galezi nauki 
i techniki. Do dziedzin, w ktorych zjawisko to ma fundamentalne znaczenie, zaliczyc 
mozna przemysl wydobywczy gazu ziemnego i ropy naftowej [D|2j|3] , hydrogeologie 
(budowazbiornikowwodnych, tarn, podziemnychujec wody) irolnictwo (odprowa- 
dzanie nadmiaru wody deszczowej i utrzymywanie odpowiedniej wilgotnosci gleby). 
Zjawisko to wykorzystuje sie tez w przemysle petrochemicznym (transport ciepla 
oraz znaczenie polimerow w procesach wydobycia surowcow) OH], przy przetwa- 
rzaniu zrodei energii (poszukiwanie materialow porowatych magazynuja_cych ener- 
gie w kolektorach sionecznych) (5), w ochronie przeciwpozarowej (transport przez 
materiaiy izolujqce od wysokich temperatur) (610, w farmacji (dostarczanie lekow 
do organizmu) (8j[9l , w ochronie srodowiska (filtracja cieczy w zbiornikach wodnych 
oraz odzyskiwanie zwiazkow zelaza z odpadow przemyslowych) HUKHKEl, diagno- 
styce medycznej (elektrody jonoselektywne uzyteczne w pomiarach stezen jonow 
w warunkach klinicznych) (13l i technologii nuklearnej (wzbogacanie uranu) fT4l . 
Transport w osrodkach porowatych lezy tez u podstaw dzialania reaktorow katali- 
tycznych, membran |14| i wiekszosci typow filtrow Stanowi tez punkt wyjscia dla 
teorii perkolacji USUI. Porowatosc i przepuszczalnosc dla ciepla, gazow i wilgoci 
to podstawowe parametry niemal wszystkich materialow budowlanych oraz mate- 
rialow codziennego uzytku (np. tekstyliow) . Jednym z ciekawszych przykladow ukla- 
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du o strukturze i cechach osrodka porowatego przepuszczalnego zarowno dla gazow 
(powietrze), jak i ph/now (krew) pluca, ktore ze wzgledu na swoja^ role w procesie 
oddychania stafy sie przedmiotem intensywnych badari doswiadczalnych i teore- 
tycznych fTTlfTBI . 

Wobec tak wielkiego znaczenia dla roznych galezi przemyslu i nauki, zagadnienie 
transportu w osrodkach porowatych jest problemem interdyscyplinarnym podejmo- 
wanym zarowno przez badaczy nauk podstawowych |T9l [201 ETJ [22] , jak i przyrodni- 
czych (21 EH, technicznych EES [26), a nawet medycznych E3EB1ESI. Mimo, ze 
w niemal niezliczonej liczbie prac zbadano juz cate spektrum roznych aspektow tego 
zjawiska i ze poswiecono mu wiele monografii (H [30l [3H EH [33l [34l [35| , wciaz kryje 
ono przed nami wiele tajemnic. 

Celem niniejszej rozprawy jest analiza dwoch szczegolowych problemow zwiq- 
zanych z transportem ph/now w osrodku porowatym. Pierwszy z nich zwiazany jest 
z zagadnieniem zwilzania materialow budowlanych. Niektore wyniki doswiadczal- 
ne i teoretyczne sugerujq, ze zjawisko to wiqze sie z tzw. dyfuzjq anomalnq. Nie- 
dawno M. Kiintz i P. Lavalee opracowali prosty model komputerowy tego procesu 
i stwierdzili, ze rzeczywiscie mamy wnim do czynieniaz dyfuzjq anomalnq. W niniej- 
szej rozprawie przedstawie alternatywne wyjasnienie wynikow uzyskanych w mode- 
lu Kiintza i Lavalee' go, nie wymagajqce odwoh/wania sie do dyfuzji anomalnej. 

Drugie zagadnienie dotyczy szczegolowej analizy numerycznej kretosci przeph/- 
wu - jednej z ciekawszych makroskopowych wielkosci fizycznych charakteryzujq- 
cych transport w osrodkach porowatych 1361137) . 

1.2. Dyfuzja normalna i anomalna 

W roku 1827 brytyjski botanik Robert Brown zauwazyt, ze pylki kwiatow znajdujqce 
sie na powierzchni cieczy poruszajq sie ruchem nieuporza_dkowanym. Zjawisko to 
pozostawalo dlugo niewyjasnione i dopiero Albert Einstein (1905) i Marian Smolu- 
chowski (1906) opisali ten ruch jako wynik mikroskopijnych zderzen pylkow z czq- 
steczkami cieczy. Odkrycie Einsteina i Smoluchowskiego ma szereg konsekwencji 
dla wspolczesnej fizyki, szczegolnie dla teorii kinetycznej oraz fizyki statystycznej. 
To dzieki nim wiadomo, ze przesuniecie sredniokwadratowe czqsteczek w transpor- 
cie dyfuzyjnym jest proporcjonalny do pierwiastka z czasu. Prawo to mozna zapisac 
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przy pomocy relacji: 

V^oc t a , (1.1) 

gdzie <r 2 > jest sredniq odleglosciq przebytq w czasie t, a wykladnik a w przypadku 
dyfuzji normalnej przyjmuje wartosc 1/2. Wielkosci makroskopowe takie jak szero- 
kosc lub pozycja profili temperatury (lub koncentracji) w transporcie dyfuzyjnym 
ciepla (lub materii) rowniez skalujq sie w czasie zgodnie z rownaniem 111.111 . tzn. opi- 
sane sajako funkcja jednej bezwymiarowej wielkosci x/vTDt. Okazuje sie jednak, ze 
istniejqprzypadki, gdy prawo to nie jest spelnione i a ^ 1/2. Mowimy wtedy o dyfuzji 
anomalnej. 

Istnieje szereg doniesien eksperymentalnych i teoretycznych na temat dyfuzji 
anomalnej pojedynczych molekut oraz anomalnej dyfuzji chemicznej (381 13911401 . 
Dotyczq one m.in. zwilzania materialow budowlanych (25j[26l|4U|42l, dyfuzji siar- 
czanu miedzi w cieczy po dejonizacji |43l|44], transportu w ukladach polimerowych 
(451 146), dyfuzji przez membrany syntetyczne S3 oraz dyfuzji w zywych ukladach 
biologicznych |48l|49] . Wpewnych warunkach anomalnq dyfuzje zaobserwowac moz- 
na rowniez w dyfuzji powierzchniowej |50|, a nawet w ruchu czasteczek w zanie- 
czyszczonej zawiesinie plazmowej (52. Fizyczne mechanizmy dyfuzji anomalnej 
w tych ukladach sq rozne i nie zawsze kompletnie poznane. Dyfuzje anomalnq poje- 
dynczych molekul tlumaczy sie miedzy innymi tym, ze odlegtosci pokonywane przez 
molekuly w kolejnych krokach czasowych lub czasy pomiedzy kolejnymi skokami 
nie spelniajq zalozeh centralnego twierdzenia granicznego i nie mog^ bye opisane 
przy pomocy rozkladu normalnego (52ll53l . 

Oprocz pojedynczych molekut, istnieje szereg ukladow, w ktorych wielkosci ma- 
kroskopowe nie spelniajq skalowania klasycznego 11.11 . a przykladem takich zjawisk 
sq procesy zwilzania (25]. Jednq z hipotez tlumaczqcych dyfuzje anomalnq w te- 
go typu ukladach jest wysunieta przez Kiintza i Lavallee'go zaleznosc wspolczynni- 
ka dyfuzji D od koncentracji dyfundujqcej substancji (54). Zasugerowali oni, ze dla 
D rosnqcego z koncentracji c wykladnik a jest wiekszy niz 1/2 i mamy do czynie- 
nia z procesami przebiegaja_cymi szybciej od klasycznej dyfuzji [ang. superdiffusion) . 
Z drugiej strony dla D malejqcego z c mamy do czynienia ze spowolnieniem tego 
procesu (44). Hipoteza ta oparta zostala na wynikach symulacji transportu pfynu 
w modelu osrodka porowatego z losowo rozmieszczonymi rozpraszaczami (551I4T1 . 
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Jednym z zagadnien poruszanych w niniejszej rozprawie bedzie krytyczna anali- 
za tezy Kiintza i Lavallee'go oraz proba podania alternatywnego wyjasnienia obser- 
wowanych przez nich efektow. 

1 .3. Kr^tosc 

Dynamika przepbywu zalezy od wielu czynnikow, poczawszy od sify zewnetrznej, po- 
przez mikrostrukture porow osrodka porowatego, predkosc przepfywu, oddziab/wa- 
nia pbynu z osrodkiem, az po fizyczny mechanizm transportu. Jednq z bardziej in- 
teresujqcych wielkosci fizycznych opisujqcych przepb/w jest kretosc - bezwymiaro- 
wa wielkosc charakteryzujqca efektywne wydluzenie drogi, wzdluz ktorej zachodzi 
transport. W klasycznym przepfywie cieczy (transportu masy) przez osrodek poro- 




a) b) 

Rysunek 1.1. Graficzne porownanie transportu przez osrodek porowaty o roznych warto- 
sciach kretosci: a) T = 1 oraz b) T> 1. 

waty kretosc T definiuje sie jako stosunek: 




(1.2) 



gdzie < A) jest srednim efektywnym wydluzeniem drogi przebywanej przez ciecz (cze- 
sto oznaczanq rowniez jako L e ), a L jest liniowym rozmiarem ukladu w kierunku 
gradientu cisnienia lub sib/ zewnetrznej, co implikuje T ^ 1 (rysunek II. 1L Kretosc 
jest wielkosciq o tyle interesuja_cq, ze zawiera w sobie dwojakq informacje: zarowno 
o strukturze osrodka, jak i o charakterze transportu. 

Juz 46 lat temu w swojej pracy opublikowanej w Nature Lorenz zauwazyt, ze kre- 
tosc w kontekscie osrodkow porowatych interpretowana jest na rozne, niekiedy wza- 
jemnie sprzeczne sposoby, co dotyczy szczegolnie jej relacji do porowatosci (561 . 
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Obecnie oprocz kretosci hydrodynamicznej w literaturze znalezc mozna odniesie- 
nia do kretosci dyfuzyjnej (57], elektrycznej (33 Qj)], a nawet kretosci w transporcie 
fali dzwiekowej [19]. Co wiecej, wyznaczane w rozny sposob zaleznosci kretosci od 
porowatosci T{<p), bardzo wazne z punktu widzenia zastosowan, rozniq sie w zalez- 
nosci od przyjetego modelu transportu oraz przyjetej metodyki badari. 

W zakresie niskich liczb Reynoldsa, gdy brak jest efektow inercyjnych (2D), prze- 
pfyw ptynu Newtonowskiego przez osrodek porowaty podlega prawu Darcy'ego: 

q=--VP, (1.3) 

gdzie q jest strumieniem pfynu, k przepuszczalnosciq osrodka, [L lepkosciq dyna- 
micznq pfynu, a VP gradientem cisnienia, pod wptywem ktorego nastepuje prze- 
ptyw Podstawowym problemem zarowno teoretycznym, jak i doswiadczalnym, jest 
zaleznosc przepuszczalnosci osrodka k od jego porowatosci (p (581 EH). Najbardziej 
znana_ zaleznosciq. tego typu jest wyprowadzona na gruncie teorii kapilarnej relacja 
Kozeny'ego ffl: 

k=c ^, (1.4) 

gdzie c jest stalq Kozeny'ego uwazanq za ceche geometrycznq ukladu, a S jest polem 
powierzchni mikroporow osrodka porowatego w jednostce objetosci (ang. specific 
surface area) . 

Historycznie, pojecie kretosci w hydro dynamice zostalo wprowadzone przez Car- 
mana w 1937 roku (60], jako poprawka empiryczna do kapilarnej teorii Kozeny'ego. 
Dzieki niej relacja miedzy przepuszczalnosci^ k a porowatosciq (p przyjela nastepu- 
ja_cqpostac: 

k = c Qj2$i> 0-.5) 
gdzie wystepujqca w tym rownaniu kretosc T, poczqtkowo wykorzystywana bardziej 
jako dodatkowy parametr teorii, okazata sie bye mierzalnq wielkosciq fizyczna,, in- 
terpretowanajako wydluzenie drogi w transporcie przez osrodek porowaty Kretosc 
nie jest jednak jedynie cechq geometrycznq osrodka, a jej wartosc zalezy rowniez 
od mechanizmu fizycznego transportu (36) . Co wiecej, jak wspomniatem wyzej, roz- 
roznia sie kretosc dyfuzyjnq, elektrycznq oraz hydro dynamicznq (SUEDES]. Kretosc 
dyfuzyjnq osrodka porowatego definiuje sie wzorem: 

D 

Td = ~, (1-6) 



gdzie to wspolczynnik dyfuzji mierzony w osrodku porowatym, a D to wspol- 
czynnik dyfuzji samego pfynu bez obecnosci przeszkod i kanalow. Ze wzgledu na 
wptyw mikrostruktury osrodka na dynamike procesow dyfuzyjnych, D ^ D^,, czyli 
I'd ^ 1 (63] . Analogicznie wprowadza sie krqtosc elektrycznq T e jako: 

T e = (p—, (1.7) 
a 

gdzie (f> to porowatosc osrodka, a a<p oraz a opornosciami wlasciwymi osrod- 
ka porowatego zanurzonego w elektrolicie oraz opornoscia_ samego elektrolitu. Ze 
wzgledu na spowolnienie procesow transportu przez obecnosc przeszkod w osrod- 
ku ^ a, czyli T e 3= 1. 

Jednym z glownych kierunkow badan nad kretosciq sq proby skorelowania jej 
z tatwiej mierzalnq doswiadczalnie porowatosciq osrodka w mozliwie jak najogol- 
niejszej formie. Wtym celu, na przestrzeni lat, przeprowadzonych zostato wiele ana- 
liz teoretycznych EUED , badan doswiadczalnych (65l[66) , a w ostatnich latach coraz 
czesciej rowniez obliczen numerycznych [671 ED, ktorych zadaniem byto wyznacze- 
nie relacji T(<p) dla roznych ukladow fizycznych. 

Najbardziej znane i charakterystyczne relacje wia_zqee T z to: 

1X0) = <p~ p , (1.8a) 

T{(p) = l-pln0, (1.8b) 

7X0) = l + p(l-0), (1.8c) 

7(0) = [l + p(l-0)] 2 , (1.8d) 
(1-0) 

7(0) = 1 + p - J (1.8e) 

(0-0 c ) m 

gdzie /? [oraz m w rownaniu lll.8el ll wolnymi parametrami uzywanymi w proce- 
durze dopasowania teorii do eksperymentu. Pierwsza z tych relacji zostala zapro- 
ponowana na podstawie badan nad przewodnictwem elektrycznym przez Archiego 
w 1942 roku |69| i jest czesto stosowana rowniez w innych kontekstach, np. do opisu 
przepfywu cieczy EH [70). Drugie z rownan zostalo wyprowadzone na gruncie teo- 
retycznych rozwazafi nad transportem dyfuzyjnym w ukladach swobodnie pokry- 
wajqcych sie sfer {p = 1/2) ED [72) lub cylindrow {p = 1 lub p = 3/2) |73|. Podob- 
na relacja zostala rowniez wyznaczona empirycznie (z parametrami p « 0.86 oraz 
p « 1.66) poprzez najlepsze dopasowanie dla kretosci hydraulicznej w eksperymen- 
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cie przepfywu przez uklady wymieszanych przeszkod z roznym stosunkiem szeroko- 
sci do grubosci (65] oraz w pomiarach kretosci elektrycznej w zawiesinie szklanych 
kul (66] . Rownanie fll .8cD jest zaleznosciq empirycznq dla przepuszczalnosci dna na- 
turalnych zbiornikow wodnych (74). Zaleznosc liniowa Jl .8cD byta rowniez wypro- 
wadzona na gruncie badan numerycznych nad przepfywem z uzyciem modelu gazu 
sieciowego {p = 0.8) |67] oraz teoretycznie dla modelu rozpraszania fali dzwiekowej 
w osrodku wypelnionym cieczq (p = 1) (6UCG). Rownanie lll.8dl l uzyskano w modelu 
kretosci dyfuzyjnej w osadach dna morskiego, dla ktorego oszacowana zostala wiel- 
kosc p « 1.1 [75]. Z kolei zaproponowana przez Koponena w 1997 roku relacja ll.8el l. 
zawiera dodatkowy wolny parametr m i zostala wyznaczona na gruncie obliczen w 
modelu zblizonym do tego, jakiego uzyjemy w analizie przeprowadzonej w niniejszej 
pracy [p = 0.65 i m = 0.19) (68). 

2.2 

2 
1.8 
1.6 
1.4 
1.2 

1 

0.2 0.4 0.6 0.8 1 

Rysunek 1.2. Porownanie zaleznosci T[(p) zebranych z kilku roznych zrodel dla roznych 
technik pomiarowych / modeli teoretycznych: a) (7U] b) (7T) c) (B7) d) [55] e) 
(64]. 

Powyzsze wzory nie sq spojne mifdzy sobq f rysunek II .2ft . Roznice wynikajq 
z takich czynnikow, jak budowa osrodka porowatego (rozny prog perkolacji), meto- 
da badawcza uzyta do wyznaczenia kretosci, rozmiar fizyczny probek przyjetych do 
analizy (efekty skonczonego rozmiaru), czy tez przyjety do badan mechanizm trans- 
portu. 
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1.4. Cele i struktura rozprawy 



Celem niniejszej rozprawy jest analiza opisanych powyzej szczegotowych proble- 
mow zwiazanych z transportem prynow w osrodku porowatym. 

Praca sktada sie z 6 rozdzialow. W rozdzialachEl i[3l opisane zostana. narzedzia 
numeryczne uzyte w trakcie badafi - automat komorkowy gazu sieciowego oraz mo- 
del gazu sieciowego Boltzmanna. Gf6wna_ czesc rozprawy stanowia_ rozdziatylH oraz 
El, gdzie przedstawione zostanq oryginalne wyniki. W rozdzialelU przedstawie ana- 
lize problemu transportu w modelu Kiinza i Lavalee'go oraz alternatywne wyjasnie- 
nie uzyskanych wynikow, nie wymagaja_ce odwoiywania sie do dyfuzji anomalnej. 
W rozdziale O przeprowadze szczegolowa^ analize numeryczna. kretosci przepb/wu 
i jej korelacji z porowatoscia, osrodka oraz z innymi wielkosciami makroskopowymi, 
ktore go charakteryzujq. Rozdziat[6l zawiera podsumowanie uzyskanych wynikow. 
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Rozdziaf 

2 



Model gazu sieciowego FHP 



W tym rozdziale przedstawie teoretyczne podstawy oraz podstawowe zasady dzia- 
tania wprowadzonego w 1986 roku przez U. Frish'a, B. Hasslacher'a i Y. Pomeaou 
automatu komorkowego FHP (55) . 

2.1. Wst^p 

Jednym z pierwszych modeli transportu ptynow, ktory nie wymaga jawnego rozwiq- 
zywania rownan rozniczkowych hydrodynamiki, jest opracowany przez Hardy'ego, 
de Pazzisa i Pomeau dwuwymiarowy model gazu sieciowego HPP (TEllZZlIZE). Jest to 
automat komorkowy o deterministycznych regulach kolizji opisujqcy dynamike dys- 
kretnych cz^steczek rozlozonych w weztach regularnej sieci kwadratowej. Ze wzgle- 
du na symetrie sieci i tylko dwa wyroznione kierunki predkosci czqsteczek na sieci, 
hydrodynamika modelu HPP jest jednak niefizyczna, gdyz otrzymywane rozwiqza- 
nia nie sq izotropowe. Co wif cej, w skali makroskopowej nie sq one niezmiennicze 
wzgledem transformacji Galileusza. W roku 1986 Frish, Hasslacher i Pomeau poka- 
zali, ze podobny automat komorkowy (nazwany FHP od inicjalow autorow), ale zde- 
finiowany na sieci trojkataej, reprezentuje w skali makroskopowej izotropowe row- 
nania przeprywu Naviera-Stokesa w zakresie niskich liczb Macha (55) . Od tego czasu 
model FHP byl przedmiotem zainteresowania roznych dziedzin nauki i techniki (79) . 
Jedna^ z bardziej interesuj^cych cech modelu jest wzgledna prostota implementacji 
skomplikowanych (w sensie geometrycznym) warunkow brzegowych, a co za tym 
idzie - mozliwosc uzycia tej metody do badan transportu w osrodkach o strukturze 
porowatej. 
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2.2. Definicja modelu 



Istnieje kilka war ian tow modelu FHP, oznaczanych skrotami FHP1, FHP2 etc. Kazdy 
z nich zdefiniowany jest na sieci trojkataej, gdzie kazdy z N x x N y weztow posiada 
szesciu sqsiadow numerowanych od 1 do 6. W kazdym wezle sieci moze sie znaj- 
dowac ... 6 cza_steczek o roznych kierunkach pedu. W zaleznosci od wariantu mo- 
del moze dopuszczac dodatkowq czasteczke spoczywaja_ca_ o pedzie rownym zeru, 
wtedy maksymalna liczba czqsteczek przypadajqca na wezel wynosi 7. W zadnym 
wezle nie mogq znajdowac sie dwie czasteczki poruszajqce sie w tym samym kierun- 
ku. W przypadkach kilku cza_steczek zmierzajqcych do tego samego wezla problem 
kolizji miedzy nimi rozwiazywany jest przez ich proste przekonfigurowanie. W za- 
leznosci od wariantu modelu ilosc i rodzaj konfiguracji, dla ktorych definiowana jest 
kolizja, rozni sie miedzy soba.- Na przyklad, najbardziej popularny wariant modelu, 
FHP3, zawiera 70 kolizji dwu- i trzy- cza_steczkowych (80j[8T|. W zaleznosci od rodza- 
ju kolizji inna jest lepkosc modelowego plynu. W niniejszej rozprawie do symulacji 
zjawisk propagacji frontu koncentracji uzyta zostanie wersja FHP5 o uproszczonych 
regulach kolizji przedstawionych na rysunku lZTl 




Rysunek2.1. Reguly kolizji dla modelu FHP5. Symbol (#) oznacza czqstke spoczywaja^, 
strzalki oznaczaj^ czasteczki poruszaja_ce sie w danym kierunku. Kompletny 
zestaw kolizji zawiera dodatkowo konfiguracje obrocone o wielokrotnosci kata 
n/3. 

2.3. Warunki brzegowe i rozpraszacze 

Uwzglednienie warunkow brzegowych w modelu gazu sieciowego nalezy do zadan 
wzglednie prostych. Wezel sieci dostepny dla cza_stek moze zostac oznaczony jako 
wezel rozpraszajqcy Implementacja warunku brzegowego bez poslizgu {ang. no slip) 
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sprowadza sie do zmiany zwrotow pedow wszystkich czasteczek znajdujqcych sie 
w wezlach rozpraszaj^cych (wiecej informacji o implementacji warunkow odbicia 
od rozpraszaczy znajduje sie w podrozdziale l4. LID . 



2.4. Algorytm 

Algorytm dzialania modelu FHP sktada sie z dwoch krokow: 1) translacji i 2) kolizji. 
Cza_steczki rozmieszczone na sieci posiadajq okreslony ped. Zgodnie z jego kierun- 
kiem w kroku translacji poruszaj^ sie one po prostych od wezta, w ktorym sie ak- 
tualnie znajdujq, do wezta sa_siedniego. Po przesunieciu czasteczek, w kroku kolizji 
nastepuje odbicie czasteczek od nieruchomych brzegow i rozpraszaczy oraz rozwia.- 
zanie kolizji miedzyczasteczkowych. To ostatnie wykonuje sie najczesciej poprzez 
umieszczenie konfiguracji „po" zderzeniu w specjalnej tablicy. 



2.5. Wielkosci makroskopowe 

Informacja, jakq otrzymujemy bezposrednio z automatu komorkowego FHP, to chwi- 
lowe konfiguracje czasteczek (ich polozenie w wezlach sieci) oraz ich pedy. W ce- 
lu wyznaczenia pol makroskopowych gestosci i predkosci w catym obszarze, nalezy 
przeprowadzic procedure usredniania. Zdefiniujmy (N a ) , jako ilosc czasteczek znaj- 
dujqcych sie w wezle i i poruszajqcych sie w kierunku a (82). Wtedy gestosc czaste- 
czek q i ped makroskopowy gu zapisac mozemy w postaci sum po wezlach i kierun- 
kach: 

gu=^=^-J^(N a c a )i, (2.2) 
V 3 M ia 

gdzie Cj jest wektorem elementarnym sieci, M jest liczba. wezlow, a czynnik -4= re- 
prezentuje ilosc wezlow na jednostke powierzchni dla sieci trojkataej o stalej sieci 
c = 1. W dalszej czesci pracy, a szczegolnie w rozdzialelH zamiast gestosci na jed- 
nostke powierzchni czesto poslugiwac sie bedziemy pojeciem koncentracji czqstek, 
ktorq definiujemyjako p = gV3/2. 
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Rozdziaf 

3 



Gaz sieciowy Boltzmanna 



W tym rozdziale przedstawie model gazu sieciowego Boltzmanna oraz omowie jego 
zastosowanie do symulacji zjawisk transportu w osrodkach porowatych. Przedstawie 
podstawy teoretyczne rownania transportu Boltzmanna z przyblizeniem jednorela- 
ksacyjnym czlonu kolizji. Omowie jeden ze sposobow rozwiazania tego rownania - 
model gazu sieciowego Boltzmanna. W celu weryfikacji opracowanego kodu nume- 
rycznego zostanie on uzyty do rozwiajzania zagadnienia przeprywu przez prostokaj:- 
ny kanal, a wyniki zostana^ porownane z rozwiazaniem analitycznym. 

3.1. Wst^p 

Automat komorkowy FHP, pomimo wielu interesujqcych wlasciwosci w kontekscie 
modelowania przeptywow przez osrodki porowate, posiada rowniez wady ktore znacz- 
nie ograniczajq obszar jego zastosowan praktycznych. Do glownych jego niedosko- 
nalosci zaliczyc mozna szum statystyczny pojawiajqcy sie w usrednianych wielko- 
sciach, bedqcy naturalnqkonsekwencj^jego mikroskopowego charakteru (83). Spo- 
sobem na obejscie tego problemu jest zwykle zwiekszenie rozmiarow obszarow, po 
ktorych usredniane sq wielkosci makroskopowe. Problem ten ma szczegolne zna- 
czenie w uktadach takich, jak osrodki porowate, gdzie mamy do czynienia z dwiema 
skalami przestrzennymi: 1) skalq mikroporow, w ktorej zachodzi transport i rozwiq- 
zywane sq rownania transportu oraz 2) skalq makroskopowq, w ktorej zaniedbuje 
sie detale i wprowadza prawa i wlasciwosci usrednione po objetosci osrodka. Juz sa- 
mo uwzglednienie obu z nich wymaga uzycia duzej siatki obliczeniowej, stqd po- 
trzeba dodatkowego jej zwiekszenia ze wzgledu na problem zaszumienia wynikow 
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w zasadzie eliminuje metode FHP z praktycznego uzycia do badania transportu hy- 
drodynamicznego w ukladach porowatych. W roku 1988 G. McNamara i G. Zanet- 
ti wprowadzili pojecie gazu sieciowego Boltzmanna (LBM) jako metody pozwalajq- 
cej wyeliminowac szum statystyczny z wynikow otrzymywanych metodami gazow 
sieciowych |83l[84]. Historycznie metoda ta wywodzi sie z klasycznych automatow 
komorkowych, w ktorych dyskretna zmienna [N a )i, oznaczajqca liczbe czasteczek 
znajdujqcych sie w wezle i i poruszajqcych sie w kierunku a (82], zostala zastajpio- 
na cia^glq zmienna. (/«)/ przyjmuja_ca_ dowolne wartosci z przedzialu [0; 1]. Wielkosc 
(f a )i oznacza w tym kontekscie prawdopodobienstwo znalezienia cza_steczki znaj- 
duja_cej sie w wezle i i poruszajqcej sie w kierunku a. Obecnie podstawe teoretyczna. 
dzialania modelu LBM stanowi teoria kinetyczna Boltzmanna |85|. 

3.2. Teoria kinetyczna 

3.2.1. Funkcja rozkfadu 

Kompletny stan uktadu mikroskopowego mozna opisac poprzez podanie zestawu 
polozen x oraz pedow p wszystkich jego cza_stek. Teoretycznie sledzenie ewolucji 
ukladu mozna zrealizowac poprzez sledzenie trajektorii kazdej z nich. Ze wzgledu 
na ich ogromnq ilosc, przy rozpatrywaniu duzych ukladow podejscie takie jest jed- 
nak niepraktyczne. Mozemy jednak przejsc z przestrzeni fazowej, w ktorej znajduja_ 
sie wszystkie czastki (atomy) reprezentujqce badany uklad, do przestrzeni konfigura- 
cyjnej polozen i pedow jednej czastki (86) . Wprowadza sie w tym celu pojecie funkcji 
rozkladu zdefiniowanej tak, by wyrazenie: 

f(x,v,t)d 3 xd 3 v, (3.1) 

opisywalo liczbe cza^steczek w skonczonym elemencie przestrzeni konfiguracyjnej 
polozen i pedow o objetosci d 3 x xd 3 v.Na rysunku lBTTI zaznaczony zostat punkt w prze- 
strzeni konfiguracyjnej wraz z otoczeniem. Wartosc funkcji rozkladu w tym punk- 
cie, /(x,v, it), mowi nam o tym, ile (srednio) cza.stek o pedzie v i polozeniu x znaj- 
dziemy w ukladzie w chwili t. Wprowadzenie funkcji rozkladu pozwolilo ograniczyc 
problemy zwiqzane z dyskretnym charakterem modelu gazu sieciowego - przede 
wszystkim wyeliminowalo potrzebe czasochtonnych usrednien wielkosci makrosko- 
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A 



d 3 x 



Jd 3 v 



Rysunek 3.1. Przestrzen konfiguracyjna poiozen i pedow z zaznaczonym jednym punktem 
wraz z jego otoczeniem o objetosci rf 3 x x rf 3 y. 

powych. Mozna wrecz powiedziec, ze proces usredniania wynikow jest tu realizowa- 
ny juz na etapie opisu ukladu przez funkcje rozkladu. 



3.2.2. Rownanie transportu 

W zagadnieniach, ktorymi bedziemy sie zajmowac, istotna jest dynamika ukladu 
oraz proces jego dochodzenia do stanu rownowagi. Dynamike ukladu w przestrzeni 
pedow i poiozen opisuje srynne rownanie transportu Boltzmanna, ktore w postaci 
cia_glej przyjmuje postac (87l : 



— + v-V r + 
[at 



^).V v )/,r,v,fl = (f) 



(3.2) 



zderz 



gdzie l^jj jest czlonem odpowiadajqcym za opis kolizji miedzycza.steczkowych, 
F - sila_ zewnetrznq dzialaja_ca_ na cza_steczki ukladu, m - masa_ cza_steczek, a f{r, v, t) 
- funkcja. rozkladu, ktorej ewolucje badamy |86| . Rownanie to bylo przedmiotem ba- 
dan na przestrzeni wielu lat, pocza_wszy od prob szukania rozwia_zafi analitycznych 
w prostych przypadkach jednowymiarowych, po badania numeryczne w dwoch 
i trzech wymiarach przestrzennych. Co ciekawe, historycznie rownanie 13.21 nie sta- 
nowilo punktu wyjscia do metody gazu sieciowego Boltzmanna. Wyprowadzenie ta- 
kie zostato przedstawione dopiero wroku 1997 (87l . 



3.2.3. Przyblizenie BGK 

Dla wyrazenia odpowiadaja.cego kolizjom w rownaniu ( 13 .2D stosuje sie najczesciej 
przyblizenie jedno-relaksacyjne opracowane przez Bhatnagara, Grossa i Krooka 
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w 1954 roku (BGK) 

(df\ f~f eq 

-77 «- — , (3.3) 

V^i /zderz T 

gdzie t to sredni czas relaksacji pojedynczego zderzenia. Funkcja rozkladu / eq opi- 
suje stan ukladu w rownowadze termodynamicznej i wyznacza sie j^ zwykle z roz- 
kladu Maxwella-Boltzmanna: 

v 2 1 

gdzie 1cb jest statq Boltzmanna, T - temperatura_ ukladu, a v - predkosci^ czqste- 
czek. Model BGK zdobyt uznanie ze wzgledu na prostote i intuicyjne potraktowanie 
zderzefi jako procesu lokalnego dochodzenia ukladu do stanu rownowagi ze stalym 
czasem relaksacji. Istnieje osobna galaz badafi nad wielorelaksacyjnymi modelami 
zderzefi w gazie sieciowym Boltzmanna - MRT {ang. multiple relaxation time) - kto- 
re odznaczajq sie wieksza^ dokladnosciq w szerszym zakresie numerycznych lepko- 
sci plynu. W literaturze znalezc mozna wiele porownan modeli BGK i MRT (8911901 . 
W niniejszej rozprawie ograniczamy sie do stosowania najpopularniejszego przybli- 
zenia BGK, ktore doskonale radzi sobie z modelowaniem zjawisk dla interesujacych 
nas parametrow 



3.3. Model D2Q9 

Rownanie 113.21 1 przybrakusilyzewnetrznej mozna zapisacwpostacidyskretnej (przyj- 
mujqc5?= 1) |HJ[83: 



//(r+e/,f+l) = / ; (r, t) + 



\dt) 



(3.5) 

zderz 



Wektory e, sq dyskretnymi wektorami predkosci rozpinajqcymi siec o zadanej syme- 
trii. Najczesciej uzywanym modelem jest wersja dwuwymiarowa o 9 mozliwych wek- 
torache/: (0,0), (1,0), (0,1), (-1,0), (0,-1), (1,1), (-1,1), (-1,-1) oraz (1,-1). Ta wer- 
sja modelu funkcjonuje w literaturze pod nazwq D2Q9, gdzie D = 2 oznacza wymiar 
przestrzenny, a Q = 9 ilosc dozwolonych wektorow predkosci. W literaturze funk- 
cjonuje szereg modeli zdefiniowanych w jednym, dwoch i trzech wymiarach prze- 
strzennych z roznq liczbq dozwolonych kierunkow wektorow predkosci 1871 . 
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3.3.1. Funkcja rozkfadu i wielkosci makroskopowe 

Znajomosc funkcji rozkladu dla badanego ukladu umozliwia wyznaczenie hydrody- 
namicznych wielkosci makroskopowych. Wielkosci opisujqce hydrodynamike bada- 
nego ukladu wyrazic mozna przez kolejne momenty funkcji rozkladu. Korzystaj^c 
z jej dyskretnej postaci i zakladajqc, ze model zdefiniowano z 9 dozwolonymi wek- 
torami predkosci, wielkosci makroskopowe takie jak gestosc g oraz predkosc lokalnq 
u mozemy wyznaczyc poprzez sumy 1921 : 

9 

i=l 
9 

gu=Y J fi^u (3.7) 

gdzie fi jest dyskretnq funkcji rozkladu, e/ sq wektorami sieci, a zakresy sumowa- 
nia w powyzszych wzorach zaleza^ od przyjetej symetrii sieci oraz liczby wymiarow 
w konkretnej wersji modelu. 

Cisnienie w modelu LBM wyznacza sie z rownania stanu: 

p = c 2 s q, (3.8) 

gdzie stala c s = l/\/3, jest predkosciq dzwieku dla omawianego modelu |93|. 



3.3.2. Model LBM w przyblizeniu BGK 

Rownanie 1 13.31 1 wraz z 113.51 1 wyrazajq ostatecznq postac dyskretnq rownania na ewo- 
lucje funkcji rozkladu w modelu gazu sieciowego Boltzmanna w przyblizeniu BGK: 

fi-f- q 



/i(r + ej,f+l) = /i(r, t) 



(3.9) 



Postac rownowagowq funkcji rozkladu otrzymac mozemy w zakresie niskich predko- 
sci, rozwijajqc rozklad Maxwella-Boltzmanna w szereg z dokladnosciq. do wyrazow 
drugiego rzedu |93|: 



9 3 
1 + 3e/ u + -(e; -u) 2 - -u 2 



(3.10) 



gdzie co = 4/9, ^1,3,5,7 = 1/9, a (o 2 ,4,6,8 = 1/36 EE |94] . 
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3.3.3. Blad dyskretyzacji 

Podstawowq kwestiq jest wielkosc btedu, jaki popelniamy rozwiazujac rownania trans- 
portu metoda. LBM, wzgledem dokladnego rozwiazania rownari Naviera-Stokesa. 
Okazuje sie E4j[84], ze rozwiazanie modelu opisanego rownaniami (13.91) w przybli- 
zeniu jednorelaksacyjnym ma dokladnosc drugiego rzedu w czasie i przestrzeni jesli 
tylko bta_d dyskretyzacji wl^czy sie w lepkosc kinematyczna. cieczy. W zwiazku z tym, 
w modelu D2Q9 lepkosc kinematyczna. definiuje sie wzorem: 

2t-1 

v=— — , (3.11) 
b 

gdzie t jest charakterystycznym czasem relaksacji w wyrazie kolizji. 

3.3.4. Warunki brzegowe 

Poprawna implementacja skomplikowanych geometrycznie warunkow brzegowych 
to jeden z najtrudniejszych problemow w obliczeniach numerycznych. Metoda LBM 
zdobyla uznanie z powodu zgodnego z intuicja_ i bardzo prostego w implementacji 
warunku na ewolucje funkcji rozkladu przy sztywnej sciance w przepbywie bez po- 
slizgu [ang. no slip boundary condition) (94). Okazuje sie jednak, ze stosowanie te- 
go prostego warunku na odbicie funkcji rozkladu przy sciankach powoduje lokalna. 
strate dokladnosci rozwia_zania i redukcje dokladnosci do pierwszego rzedu. Dzieki 
zastosowaniu statego czasu relaksacji t = 1 uzytego w niniejszej pracy i zastosowaniu 
ulepszonej procedury odbicia {ang. halfway bounce back) uzyskalismy dokladnosc 
drugiego rzedu w czasie i przestrzeni |84 | . 

3.3.5. Algorytm 

Typowy algorytm metody gazu sieciowego Boltzmanna sprowadza sie do wykonania 
dwoch krokow: 

1) kolizji, czyli uwzglednienia czlonu kolizji zgodnie z rownaniem: 

fdr,t)-f°Hu,Q) 
fi(r, t) = fi(r, t) ! , (3.12) 

T 

gdzie fi(r, t) to pewne zmienne pomocnicze, 
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2) propagacji, czyli transportu funkcji rozktadu zgodnie z rownaniem Boltzman- 
naEHll: 

fi(r+e i ,t+l) = f i (r,t). (3.13) 

Algorytm symulacji sprowadza sie wiec do rozdzielenia zagadnienia transportu na 
dwa oddzielne procesy - transportu i lokalnych kolizji cza_stek, czyli tak, jak w przy- 
padku omawianego wczesniej modelu FHP. 



3.4. Weryfikacja kodu numerycznego 

W celu weryfikacji implementacji metody LBM stworzonej na potrzeby niniejszej 
rozprawy, rozwiazany zostal problem przepb/wu przez prosty kanat w obecnosci si- 
ly wymuszaja_cej (lub rownowaznej roznicy cisnien na wlocie i wylocie z kanalu). Ze 
wzgledu na istniejqce rozwiazanie analityczne problem ten, znany wliteraturze jako 
przeplyw Pouseuille'a [ang. Poiseuille flow), jest czesto uzywany do weryfikacji po- 
prawnosci implementacji numerycznej. Na rysunku l3~2l przedstawiony zostal sche- 
mat omawianego problemu. Obszar prostoka_tny ograniczony punktami xq, X\, yo 



yi 



yo 




Pout 



H > 

X! X 

Rysunek 3.2. Schemat problemu przepfywu przez zamkniety kanal {ang. Poiseuille flow). 

oraz yi jest wypelniony prynem. Scianki poziome w polozeniach y = yo oraz y = yi 
sq sztywne (brak poslizgu, t.j. predkosc v x = na sciankach). Ze wzgledu na roznice 
cisnien pomiedzy przekrojami kanalu dla x = Xq oraz x = x\ w obszarze kanalu na- 
stepuje przeplyw cieczy. Obliczenia analityczne i rozwiazanie uproszczonego row- 
nania przeplywu Naviera-Stokesa prowadzi do wniosku, ze profil predkosci w takim 
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-+- 

Xo 



kanale ma charakter paraboliczny i przy zalozeniu, ze yn = i y\ = h spelnia rowna- 
nie: 

v(x,y) = {v x , v y ) 



( 1 dp \ 

o,-—-fy{y-h)\, (3.14) 
^ 2vp dx ) 



gdzie jest gradientem cisnienia miedzy wlotem, a wylotem z kanaiu, v - lepkosciq 



2vq dx 
lotem, a 

kinematycznq, g - gestosciq ptynu, ah- szerokoscia. kanaiu 



Warunki brzegowe i poczqtkowe 

Zgodnie z oznaczeniami z rysunku 13.21 zaktadamy warunki pocza_tkowe predkosci 
oraz cisnienia w postaci: 

v x {x,y, t = 0) = i/ y (x,y, £ = 0) = 0, 

y (3.15) 
p{x,y, t = 0) = p , 

w cafym obszarze kanaiu, gdzie po = (pi n - p ut)/2 {ang. in - wlot, ang. out- wylot 
z kanaiu) . 

Jako warunki brzegowe cisnienia na wlocie [x = x ) i wylocie z kanaiu (x = X!) 
w dowolnej chwili przyjete zostaty stale cisnienia p m i p out , gdzie p m > p out . Komplet 
warunkow brzegowych, wlaczajqc te nakladane na pole predkosci, wyrazic mozna 
nastepujqco: 

v x (x,y Q ,t) = v x {x,yi,t) = v y {x,y ,t) = v y (x,yi,t) = 0, 

p{x Q ,y,t) = p m , t^O, (3.16) 
p{x x ,y, t) = p ouU t^O. 

Porownanie z rozwiazaniem analitycznym 

W zagadnieniu testowym, pomiedzy wlotem, a wylotem z rury, przyjety zostal gra- 
dient cisnienia Ap = p in - p out = 0.1 mu ts~ 2 {ang. mass unit per time step squared) 
(96) . Problem zostal rozwiazany na sieci prostokataej o wymiarach 32 x 16 dla lepko- 
sci kinematycznej v = 1.0 (co daje t = 0.4). Przyjmujqc yo = 0, y\ = 1, xo = i Xx = 1, 
rozwiazanie analityczne pol predkosci i cisnienia dla omawianego problemu mozna 
zapisac wwygodnej postaci: 

p{x,y, t) = p m -a-x, 

v x {x,y,t) = b-y{\-y), (3.17) 
v y {x,y, t) = 0, 
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Rysunek3.3. Porownanie profilu cisnienia dla problemu przepfywu przez prosty kanal 
z rozwiazaniem analitycznym 13.17t dla a = 0.1. 
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Rysunek3.4. Porownanie profilu predkosci dla problemu przeplywu przez prosty kanal 
z rozwiazaniem analitycznym Q.17) dla b = 1.77. 



gdzie zmienne a oraz b uzyte do dopasowania wzorow do danych pomiarowych 
(97). Na rysunkach 13.31 oraz 13.41 przedstawione zostary wyniki odpowiednio dla ci- 
snienia oraz predkosci w kierunku x odczytanej wzdluz przekroju przez kanal. Pa- 
raboliczny profil predkosci na przekroju kanatu oraz liniowa zaleznosc cisnienia od 
odleglosci do wlotu rury uzyskane przy pomocy modelu LBM odpowiadaja, dosko- 
nale przedstawionym rozwiqzaniom analitycznym. 
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Rozdziat 

4 

Transport w modelu 
Kuntza-Lavallee 'go 

W tym rozdziale przedstawie krytycznq analize wystepowania dyfuzji anomalnej 
w modelu Kiintza-Lavallee'go (KL). Uwzgledniona i zbadana zostanie rola i znacze- 
nie dryfu hydro dynamicznego w powstawaniu obserwowanych zjawisk „anomal- 
nych". 

Rozdziat ten opiera sie na wynikach opublikowanych w pracy: 
M. Matyka, Z. Koza, 

Spreading of a density front in the Kuntz-Lavallee model of porous media, 
J. Phys. D: Appl. Phys. 40, 4078-4083 (2007). 

4.1. Opis zagadnienia 

Jednq z podstawowych wtasciwosci klasycznej dyfuzji pojedynczej czqsteczki jest 
skalowanie sie jej przesuniecia sredniokwadratowego {ang. root-mean-square displa- 
cement) \J (r 2 ) z czasem: y (r 2 > oc t a , gdzie a = 1/2. Istniejq jednak doniesienia 
o wielu naturalnych zjawiskach, w ktorych skalowanie z wykladnikiem 1/2 nie za- 
chodzi, a w ostatnich latach wysunieto nawet hipoteze, iz zjawisko dyfuzji anomal- 
nej wystepuje w uktadach, w ktorych wspotczynnik dyfuzji zalezny od koncentracji 
dyfundujqcej substancji |54|. Hipoteza ta zostala uzasadniona zarowno przy uzyciu 
symulacji modelem gazu sieciowego, jak i na gruncie doswiadczalnym EH |98] . Za- 
sadniczym celem dalszej czesci tego rozdzialu jest weryfikacja powyzszej hipotezy 
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powtorzenie otrzymanych wczesniej wynikow oraz proba ilosciowego wytlumacze- 
nia otrzymanych odchyleri od wykladnikowklasycznych bez odwoh/wania sie do dy- 
fuzji anomalnej. 



4.1.1. Model Kuntza-Lavallee'go 

Model Kuntza-Lavallee'go (KL) symulacji transportu ph/nu przez osrodek porowaty 
oparty zostal na standardowym modelu gazu sieciowego FHP. Ptyn reprezentowany 
jest przez jednakowe, dyskretne czasteczki rozlozone w wezlach sieci o symetrii hek- 
sagonalnej. Wkazdym wezle sieci znajdowac sie moze maksymalnie do 7 czqsteczek, 
kazda o innym kierunku predkosci (wlqezajac w to czastki spoczywajqce o pe dzie 
p = 0). Do reprezentacji mikrostruktury osrodka porowatego wykorzystane zostaly 
komorki rozpraszajqce rozlozone losowo w obszarze, w ktorym zachodzi transport 
substancji (rysunek l4.1l l, co odroznia model KL od standard owego gazu sieciowego. 
Pojedynczy krok symulacji sklada sie z dwoch czesci: 1) translacji wzdluz wektorow 



y A 



C1 



V 



c 2 



Rysunek 4.1. Model osrodka porowatego zbudowanego z rozpraszaczy (■). Linia przerywa- 
na to schematyczny profil koncentracji pocza_tkowej z uskokiem przy x = 0. 



sieci oraz 2) kroku kolizji. Pierwszy jest standardowym balistycznym ruchem czq- 
steczek wzdluz wektorow sieci. W drugim kroku, oprocz kolizji czqstek miedzy sobq 
(opisanych wczesniej w rozdzialeEJ, uwzglednione sq odbicia od rozlozonych w ob- 
szarze transportu rozpraszaczy. Odbicie to zachodzi calkowicie sprezyscie (czastecz- 
ki zachowujq swojq energie, zmieniajqc jedynie kierunek pedu). Przyklad rozwiaza- 
nia kolizji czasteczki z rozpraszaczem zostal przedstawiony na rysunku l4~Tl 
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Rysunek4.2. Przyklad kolizji cz^stek z rozpraszaczem (■) umieszczonym w wezle sieci. 

Zgodnie z zasadq zachowania energii, czqstka odbija sie od przeszkody, zmie- 
niaj^c kierunek pedu i zachowuj^c jego wartosc. 

4.1.2. Warunki brzegowe 

W celu zbadania dynamiki rozchodzenia sie frontow gestosci, na koncentracje pfynu 
nalozony jest warunek schodkowy reprezentowany rownaniem: 

c{x, t 3= 0) = C\H{-x), x < 0, 

(4.1) 

c{x, t-0) = c 2 H{x), x>0, 

gdzie H{x) jest funkcja. skokowa^ Heaviside'a, a punkt x = jest granicq oddzielajqcq 
obszary o koncentracji poczadcowej C\'\C2- Za warunek pocza_tkowy rozpatrywanego 
problemu przyjmiemy iloraz koncentracji C\/cz > 1. Dzieki istnieniu gradientu kon- 
centracji na granicy x = 0, spodziewany jest transport czasteczek gazu w kierunku 
obszaru o nizszej koncentracji. Na gornej oraz dolnej sciance ukladu zastosowane 
zostafy warunki periodyczne. Scianki pionowe ustalone zostaty jako sztywne prze- 
szkody od ktorych czqstki odbijajq sie sprezyscie. Zgodnie z rownaniem II4.1L w ob- 
szarze x < utrzymywana jest koncentracja c[x) = C\ dla wszystkich t > 0, dzieki 
czemu zapewniony jest staty dopfyw czasteczek gazu do ukladu. Dzieki uzyciu siatek 
obliczeniowych o duzej dlugosci uktad zachowuje sie jak uktad nieskonczony (brak 
widocznego wprywu prawego brzegu na otrzymywane rezultaty). 
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4.2. Wyniki 



4.2.1. Propagacja frontu koncentracji 

Przeprowadzilismy symulacje rozchodzenia sie frontu koncentracji na sieci o roz- 
miarach LXx LY = 10000 x 1000 z koncentracji rozpraszaczy rowna_ c s = 0.08. Warun- 
ki brzegowe zostaty ustalone zgodnie z rownaniem 1 14.11 1, w ktorym przyjeto C\ = 0.9 
i c 2 = 0.2. Rozktady cza_stek dla t = 1600 x 2 k , k = 0, 1,2,3 zostary przedstawione na 
rysunku l4.3i 




Rysunek4.3. Propagacja frontu gestosci dla czasow t = 1600 x 2 k , k = 0, 1,2,3. Ciemniejsze 
obszary odpowiadaj^ wiekszej koncentracji czastek. 
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Profile gestosci dla t = 8000 x 2 k , k = 0, 1, . . . , 6 zostafy pokazane na rysunku l4~4l 
w pol-logarytmicznym uktadzie wspolrzednych. Ksztalty oraz odleglosci pomiedzy 




100 1000 10000 

x 

Rysunek4.4. Profile koncentracji c[x, t) dla czasow t = 8000 x 2 k , k = 0, 1,...,6. Warunki 
brzegowe koncentracji C\ = 0.9 i C2 = 0.2. 




Rysunek4.5. Profile koncentracji c(x. t) z rysunku !4.4l w funkcji xt 112 . 

krzywymi koncentracji dla poszczegolnych pomiarow nie rozniq sie bardzo miedzy 
sobq, co pozwala wysunqc hipoteze, ze profil c{x, t) mozna opisac jako funkcje jed- 
nej zmiennej xl t a . Mozna to tatwo zweryfikowac, rysujqc profile koncentracji dla 
roznych czasow w funkcji xl t a . Jesli hipoteza o skalowaniu jest prawdziwa - profi- 
le powinny ulozyc sie jeden na drugim w pojedyncza_ krzywq. W naszym przypadku 
proba skalowania profili z wykladnikiem a = 1/2 nie jest udana, co widac wyraznie 
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na rysunku l4~5l (profile nie ukladajq sie na jednej krzywej). Analiza taka zostata prze- 
prowadzona wczesniej w pracy |54|, w ktorej pokazano, ze najlepsze dopasowanie 
uzyskuje sie dla a ~ 0.55, co zinterpretowano jako super-dyfuzje o a > 0.5. Zaloze- 
nie o super-dyfuzji powoduje jednak dose paradoksalny wynik przyspieszenia trans- 
portu poprzez dodanie rozpraszaczy. Ktintz i Lavallee w (54] thimaczyli to na dwa 
sposoby. Pierwszy to super-dyfuzyjny charakter transportu w ukladach, w ktorych 
wspolczynnik dyfuzji zalezy silnie od koncentracji substancji (z takim wlasnie ukla- 
dem mamy do czynienia w przypadku modelu FHP). Drugie wytlumaczenie mowi, 
ze anomalna dyfuzja jest tylko zjawiskiem przejsciowym i dla dostatecznie duzych 
czasow t — ► oo wykladnik a bedzie bardzo wolno zbiegal do klasycznej wartosci 1/2. 

Jednak zgodnie z analizq Botzmanna-Matano 199111001 . kazde rozwiqzanie row- 
nania dyfuzji postaci: 

D{c)—^\, (4.2) 



dt 



dx 



dx 



z warunkami brzegowymi ( 14. ID skaluje sie z wykladnikiem a = 1/2 dlakazdego t, nie- 
zaleznie od postaci D{c) 11001 . W celu weryfikacji tego twierdzenia rozwiazalem nu- 
merycznie rownanie 114.211 z uwzglednieniem dokladnej postaci funkeji D{c) dla mo- 
delu FHP5 bezposrednio z (54). Wynik zostat przedstawiony na rysunku l4!6l Widac 

100 c — 



symulacia • 
f(t)=at a5 



o 
ii 



10 



10 



100 



1000 



10000 



t 



Rysunek 4.6. Zaleznosc polozenia wartosci c = 0.5 frontu koncentracji w funkeji czasu. Roz- 
wia_zanie numeryczne rownania 14.2t ze wspolczynnikiem D(c) odczytanym 
z pracy (54) . 



wyraznie, ze otrzymane rozwi^zanie rownania dyfuzji skaluje sie klasycznie w cafym 
obszarze badanych czasow. 
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W zwiazku z powyzszymi uwagami, brak skalowania frontu koncentracji z wy- 
ktadnikiem 1/2 w modelu KL oznacza, ze rownanie klasycznej dyfuzji 114.21 1 nie moze 
skutecznie opisac tego zjawiska. Oprocz dyfuzji, w obserwowanej propagacji fron- 
tu koncentracji, musi brae udziai rowniez inny mechanizm transportu powodujqcy 
przyspieszonq migracje cza_steczek z obszaru o wyzszej koncentracji wglab osrod- 
ka porowatego. Aby zrozumiec dynamike transportu, wyznaczona zostala srednia 
predkosc (usredniona po cafym obszarze symulacji) przypadajqeana cza_steczke pfy- 
nu odpowiednio dla kierunkow x oraz y. Na rysunku |4~71 przedstawiona zostata za- 
leznosc tych wielkosci od czasu w badanym uktadzie. Widac wyraznie, ze zgodnie 
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Rysunek 4.7. Predkosc kierunku x oraz y przypadaj^ca na jednq czqsteczkg piynu od czasu. 

Wartosci podane zostaiy w jednostce lu tu~ l (stala sieci / krok symulacji). 

z oczekiwaniami, predkosc cieczy w kierunku y oscyluje wokol 0. Inaczej jest w przy- 
padku predkosci w kierunku x, ktora rosnie bardzo szybko do wartosci okolo 3 • 10 -3 
lu tu~ l (stata sieci / krok symulacji), by nastepnie bardzo powoli malec, utrzymu- 
jqc wartosc dodatniq. W poprzednich badaniach innych autorow ED wptyw nieze- 
rowego pedu cza.steczek zostal pominiety (ze wzgledu na jego pozornie niewielkie 
wartosci). W dalszej czesci tego rozdzialu bede staral sie pokazac, ze efekt ten jest 
istotny Sprobuje odpowiedziec na pytanie, czy poszukiwanym mechanizmem od- 
powiadaja_cym za anomalne zachowanie moze bye hydro dynamiczny przepfyw cie- 
czy spowodowany wysokim gradientem cisnienia wystepujqeym pomiedzy obszara- 
mi o roznych koncentracjach. 
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4.2.2. Srednia pr^dkosc frontu 

W celu sprawdzenia poprawnosci hipotezy o przeprywie hydrodynamicznym w ukla- 
dzie, przeprowadzona zostala analiza ilosciowa przestrzennego rozkladu sktadowej 
x wektora predkosci, v x {x, t), usrednionego po wszystkich komorkach w danej ko- 
lumnie. Na rysunku l4~8l przedstawione zostaly profile predkosci dla trzech roznych 
chwil czasu. Z rysunku widac wyraznie, ze dla kazdego z czasow obszar x > mo- 
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Rysunek 4.8. Skladowa x wektora predkosci - v x {x, y, t) - usrednionego po y dla trzech roz- 
nych chwil czasu t = 8000 x A k , k = 0, 1, 2. 

ze zostac podzielony na dwa charakterystyczne podobszary. W jednym z nich ped 
przeph/waja_cej substancji jest mniej wiecej stary i jest to obszar penetracji substan- 
cji w glab osrodka porowatego. To w tym obszarze uklad odpowiada na istniejqcy 
gradient cisnienia zwiekszeniem pedu. Druga czesc uktadu nie zostala spenetrowa- 
na przez wph/waja_ca_ substancje i srednie wartosci v x {x, t) oscyluja, wokol zera, co 
odpowiada warunkowi pocza_tkowemu. Z wykresu widac, ze z uph/wem czasu war- 
tosc maksymalna pedu w obszarze frontu maleje, zwieksza sie za to szerokosc tego 
obszaru. W celu ilosciowej analizy dynamiki propagacji frontu, zdefiniuje predkosc 
frontu koncentracji vr. Wogolnosci wielkosc tajestfunkcja, x i t, ale nas interesowac 
bedzie jej przyblizona wartosc w czasie t. Mozna ja. oszacowac jako srednia, wazona_ 
z waga, c[x, t) - c 2 wzdluz kierunku propagacji frontu: 

poo 

J v x {x, f)[c{x, f)-c 2 ]dx 

W) = ^—735 • (4.3) 

I [c{x,t) - c 2 ]dx 
Jo 
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Dzieki wykorzystaniu usredniania wazonego z wagami pochodzqcymi bezposred- 
nio z profili koncentracji znacznie zmniejszony zostal wpfyw szumu pochodza_cego 
z modelu gazu sieciowego w obszarach niskich koncentracji, gdzie c[x, t) ~ c 2 (rysu- 
nek !4.81 x > 3000) . Korzystajqc z definicji wyrazonej wzorem 1 14.31 1 wyznaczona zosta- 
la zaleznosc i)f(t) frysunek !4.9L Z rysunku widac wyraznie i ilosciowo, ze predkosc 
srednia frontu w chwilach poczqtkowych jest bardzo wysoka i bliska wartosci mak- 
symalnej dla modelu sieciowego Vf = 1. Dla wiekszych czasow Vf maleje jak funkcja 
wykiadnicza t~P z wyktadnikiem /3 < 1/2 bardzo siabo zalezqcym od czasu. 
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Rysunek 4.9. Predkosc Df frontu koncentracji w funkcji kroku czasowego. Dla orientacji za- 
znaczona zostala linia o nachyleniu z wykladnikiem -1/2. 

Znajomosc sredniej predkosci frontu koncentracji pozwolila na wyznaczenie sred- 
niej drogi s[t) przebytej przez front w czasie t, jako bezposredniej konsekwencji ist- 
nienia niezerowej predkosci (przepfywu): 

slt)= [ DfMdT, (4.4) 
Jo 

ktora dla duzych czasow rosnie jak t l ~^ l \ czyli szybciej niz byloby to w przypadku 
klasycznej dyfuzji z wykladnikiem 1/2. Bezposrednio z wynikow numerycznych od- 
czytac mozna pozycje Xf frontu dla ktorej c(Xf) = 0.5, ktora wynosi Xf « 3750 lu (jed- 
nostek sieci) po czasie t = 256000 tu (jednostek czasu) frysunek l4~4l . Jednoczesnie 
calka 1 14.41 1 dla tego samego t daje s{t) « 1440 lu. Oznacza to, ze okolo 38% z przesu- 
niecia x f bylo spowodowane przeplywem cieczy, a reszta dyfuzja_. Dlatego przepfyw 
hydrodynamiczny nie moze bye pominifty w analizie dynamiki ukladu. 
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4.2.3. Ruchomy ukfad odniesienia i skalowanie 

W celu oddzielenia mechanizmu transportu hydrodynamicznego od dyfuzyjnego 
przeprowadzone zostato klasyczne przejscie do ruchomego ukladu odniesienia {%', t), 
gdzie: 

x' = x-s{t). (4.5) 

Spodziewany rezultat takiej transformacji to odseparowanie lub zminimalizowanie 
efektow pochodz^cych od przeptywu z niezerowq predkosciq i obserwacja klasycz- 
nej dyfuzji z wykladnikiem 1/2. Jesli powyzsza hipoteza o klasycznej dyfuzji bylaby 
prawdziwa, profile koncentracji w ruchomym ukladzie odniesienia powinny spel- 
niac relacje skalowania: 

c{x',t) = f{(x'-x )ly/t), (4.6) 

gdzie / jest funkcjq podobienstwa, a Xo jest stal^. Korzystaj^c z wartosci pozycji fron- 
tu dla c{x' , t) = 0.5 oszacowana zostala staia x « -190, a nastepnie profile koncen- 
tracji c w funkcji (x? - Xo) lyft wyznaczono dla t = 8000 x 2 k , k = 0, 1, . . . , 6 (rysunek 
I4.10II . Na rysunku widac wyraznie, ze przejscie do ruchomego ukladu odniesienia 
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Rysunek 4.10. Profile koncentracji c w funkcji (jc 7 -Xq) l\ft dla Xo = - 190. Wstawka pokazuje 
centralnq czesc krzywych w osiach nielogarytmicznych i powiekszeniu. 

pozwolilo wyodrebnic klasycznq dyfuzje, a krzywe profili koncentracji asymptotycz- 
nie ukladajq sie na jednej krzywej. Dla czasow wiekszych niz t = 64000 przeskalowa- 
ne profile sq praktycznie nierozroznialne. 
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4.2.4. Nieciaglosc profili przy x = 



Oddzielenie zjawiska przepfywu od dyfuzji w modelu KL pozwala rowniez wyttu- 
maczyc kolejna_ z jego charakterystycznych cech - nieciaglosc profilu koncentracji 
w okolicach x = 0. Mimo ze warunki brzegowe wyrazone rownaniem (14. ID utrzymujq 
stala, koncentracje c\ w obszarze x < 0, obserwuje sie niecia_glosc koncentracji i dose 
znacznie obnizonq wartosc (duzo nizsza. niz c\) zaraz za punktem brzegowym (rysu- 
nek !4.4l dla C\ = 0.9). Niecia_glosc koncentracji Ac jest zwiazana z nieciqgloscia, pred- 
kosci sredniej v x przepb/wu na brzegu v x = dla x < i v x « Of dla x > 0. Oznaczmy 
koncentracje w punkcie x = 1 lu (stafych sieci) jako c + = C\ - Ac. Tempo transportu 
cza_stek przez piaszczyzne x = jest w przyblizeniu rowne 3ci - 3c + = 3Ac (czynnik 
3 jest liczba. mozliwych kierunkow z ktorych moze zachodzic transport w kierunku 
poziomym ze strony lewej na prawaj. Wartosc ta musi bye zrownowazona przez hy- 
drodynamiczny strumien czqstek, ktory moze bye przyblizony przez 7c + Vf (czynnik 
7 jest liczba. mozliwych kierunkow wektora predkosci). Sta_d, 3ci - 3c + « 7c + i>f, co 
pozwala oszacowac wielkosc nieciqglosci c: 

7cii)f 



Ac 



(4.7) 



3 + 7i> f 

Porownanie wyprowadzonej relacji z wynikami z symulacji zostalo przedstawione 
na rysunku Klll 
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Rysunek 4.11. Wielkosc spadku koncentracji Ac na brzegu x = w funkeji predkosci frontu 
Df, symulacja (□) i teoria (-). 
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4.3. Omowienie wynikow 



Wiele wlasciwosci modelu uzytego w tym rozdziale do badania propagacji frontow 
koncentracji zasluguje na blizsza. analize. Po pierwsze, do reprezentowania osrodka 
porowatego i redukcji efektow hydro dynamicznych w modelu uzyte zostafy rozpra- 
szacze punktowe (54). Dzialanie tych rozpraszaczy jest jednak osobliwe: sktadowa 
x predkosci zmienia sie jedynie przy kolizji z rozpraszaczem cza_stek poruszajqcych 
sie dokladnie wzdluz kierunkow rownolegfych do osi x. W omawianym modelu KL 
rozpraszacze zostah/ rozlozone losowo w obszarze osrodka porowatego. Okazuje sie 
wiec, ze taka konfiguracja nie ma duzego wptywu na dyfuzyjnosc w kierunku osi x. 
Wiadomo rowniez, ze srednia droga swobodna cza.steczek, a co za tym idzie dyfuzyj- 
nosc w oryginalnym modelu FHP (bez rozpraszaczy), dazy do oo z koncentracji da_- 
zqca. do lub 1 (54) . Biora_c pod uwage dwa wyzej wspomniane fakty latwo jest wytlu- 
maczyc zaleznosc wspolczynnika dyfuzji od koncentracji cza_steczek, przedstawionq 
wczesniej jako argument na rzecz rzekomej anomalnej dyfuzji transportowanej sub- 
stancji (54). Model KL jest oparty w glownej mierze o hydro dynamiczny model gazu 
sieciowego FHP w ktorym lokalne zasady kolizji i translacji spelniajq zasady zacho- 
wania pedu i masy i prowadzq do rownah przephywu w skali makroskopowej. Jedna. 
z cech tego typu modeli jest przeph/w od obszarow o wyzszym do obszarow o niz- 
szym cisnieniu, przy czym w modelu FHP cisnienie jest proporcjonalne do koncen- 
tracji (781 . 

Mozna oszacowac skale efektow hydrodynamicznych. Zakladajqc niskq predkosc 
frontu koncentracji, strumien czqsteczek q przeph/waja_cych przez osrodek porowaty 
jest proporcjonalny do gradientu cisnienia, co wyraza prawo Darcy'ego[H: 

q = --VP, (4.8) 

gdzie k jest stalq przepuszczalnosci (cecha osrodka porowatego), [i jest lepkoscia, dy- 
namiczny cieczy, a VP jest gradientem cisnienia miedzy wlotem a wylotem cieczy. 
Skoro (70c v x [2,jucx D~ l (przyblizenie Stokesa-Einsteina oph/wu cieczy wokot ku- 
li) oraz Poc c (80], wiec: 

v x oc DVc. (4.9) 
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Z powyzszego wzoru wynika, ze dla duzych wartosci dyfuzyjnosci D (ktore, jak wspo- 
minalem wyzej, wystepujq w modelu KL) nawet niewielki gradient koncentracji Vc 
czqstek moze indukowac znacznq predkosc unoszenia. 

Kolejnq ciekawa^ cecha^ w badanym modelu jest duza wartosc stalej Xo w rela- 
cji II4.6L Co wiecej, czas po ktorym relacja ta jest spelniona jest rowniez dose du- 
zy. Mozna podejrzewac, ze oba te efekty sq wynikiem ograniczenia z gory predko- 
sci v x ^ 1 i relacji wyrazonej rownaniem D4.9I I. Z warunkow brzegowych wiadomo 
bowiem, ze w obszarze frontu koncentracji w chwilach poczadcowych panujq wyso- 
kie wartosci Vc oraz D, natomiast v x jest ograniczone z gory przez 1 lultu. Dlatego 
prawo Darcy'ego moze bye spelnione dopiero po pewnym czasie relaksacji ukladu 
potrzebnym do redukeji gradientu Vc do wartosci rzedu 1 ID, co odpowiada warun- 
kowi na maksymalnq wartosc v x « 1 lultu. To spostrzezenie tlumaczy dlaczego w po- 
przednich doniesieniach anomalnq propagacje frontow koncentracji obserwowano 
w przypadku, gdy C\ byto bliskie 1, a rozpraszacze miary forme przedstawionq w tym 
rozdziale EH HU [98) . Warto zauwazyc, ze dose duza wartosc xq oraz duze wartosci 
czasu relaksacji ukladu sq zgodne z aktualnymi badaniami nad zwilzaniem mate- 
rialow budowlanych, gdzie krotkozasiegowe odchylenia od skalowania z czasem jak 
t 112 zostaty zaobserwowane juz wczesniej 1421 . Co wiecej, z rysunku 1431 wynika. ze 
wartosci predkosci Of dla mah/ch czasow sq duzo mniejsze niz wartosci oczekiwane 
z ekstrapolacji prawa Darcy'ego. Konsekwencja^ tego jest ujemna wartosc x i „ano- 
malny" transport profilu koncentracji w poczqtkowej fazie. Mozna sie spodziewac, 
ze w ukladach w ktorych ekstrapolacja wartosci predkosci frontu daje mniejsze war- 
tosci od tych mierzonych, wartosc xo bylaby dodatnia, a transport profili zachodzil 
z wykladnikiem a < 1/2. Podsumowujqc, z rysunku 1431 wynika, ze dla dostatecznie 
duzych czasow symulacji (t > 10 3 ) predkosc srednia frontu Vf oc l/\/7, co prowadzi 
do oszacowania charakterystycznej dtugosci hydrodynamicznej J T f =0 1/v^dr oc y/t. 
Prowadzi to do wniosku, ze skale dtugosci zjawisk dyfuzyjnych i hydro dynamicznych 
w omawianym modelu sq takie same i nie jest mozliwe rozroznienie pomiedzy nimi 
na gruncie badan nad asymptotycznym zachowaniem frontow koncentracji. 



35 



4.4. Podsumowanie 



W tym rozdziale udalo sif pokazac, ze w modelu gazu sieciowego FHP uzytego do 
analizy propagacji frontu ge stosci oprocz klasycznej dyfuzji Ficka wystepuje nieze- 
rowy ped materii, a polozenia profili gestosci skalujq sie zgodnie z prawem pote- 
gowym jak f 2 (gdzie t oznacza czas) w ruchomym uktadzie odniesienia. Podejscie to 
okazato sie pomyslem prostszym i bardziej efektywnym niz opublikowane wczesniej 
w literaturze doniesienia o dyfuzji anomalnej (HUSH [98), pozwalaja_c wytlumaczyc 
m.in. spadek koncentracji przy brzegu x = 0. Analiza przeprowadzona w tym roz- 
dziale pozwala spojrzec calosciowo na przebieg transportu w modelu KL. Nalozony 
warunek brzegowy - poczqtkowy gradient koncentracji wzdluz brzegu x = - jest 
rownowazny gradientowi cisnienia. W konsekwencji prowadzi to do powstania sit 
wymuszaja_cych przeptyw. Przepb/w ten jest blokowany przez rozlozone losowo roz- 
praszacze, sta_d mozna sie spodziewac, ze wtakim ukladzie efekty hydrodynamiczne 
beda_ asymptotycznie (w czasie) zaniedbywalne. Okazalo sie jednak, ze sposob roz- 
praszania, wyroznienie kierunkow rownolegh/ch do osi x oraz nieefektywnosc roz- 
praszaczy dla koncentracji bliskich 1 powoduja., iz w modelu wystepuja. efekty hy- 
drodynamiczne, ktore wpfywaja, na dynamike frontow koncentracji. Efekty te maja. 
istotny wph/w na czas relaksacji ukladu do stanu, w ktorym spelnione jest prawo 
Darcy'ego, a co za tym idzie - na przesuniecie x frownanie 14.61 1 . Propagacja fron- 
tu koncentracji jest wiec realizowana przy pomocy dyfuzji oraz przeph/wu cieczy 
przy czym dla dostatecznie duzych czasow oba procesy daja_ przesuniecie propor- 
cjonalne do czasu jak y/t. W konsekwencji dla dostatecznie duzych czasow profile 
koncentracji rnoga, zostac opisane jako funkcje pojedynczej zmiennej xl \ft, co moze 
prowadzic do falszywego wniosku o wylacznie dyfuzyjnym charakterze omawianego 
procesu. 
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Rozdziaf 

5 



Krptosc przeptywu 



W rozdziale tym przedstawie badania nad kretosci^ T przeptywu cieczy przez mikro- 
skopowy model osrodka porowatego. Wprowadze definicje kretosci jako bezwymia- 
rowego parametru fizycznego opisujqcego stopieri wydluzenia torow wybieranych 
przez ph/n podczas przeph/wu. Wyznacze numerycznq zaleznosc T od porowato- 
sci <p z uwzglednieniem efektow skoficzonego rozmiaru sieci. Wyniki porownam do 
zaleznosci empirycznych. Rozdzial ten opiera sie na wynikach opublikowanych w 
pracy: 

M. Matyka, A. Khalili, Z. Koza, 
Tortuosity-porosity relation in porous media flow, 
Phys. Rev. E 78, 026306 (2008). 

Dodatkowo, w podrozdziatach 15.2. 1L 15.4.51 oraz 15.4.61 przedstawione zostana. nowe 
wyniki dotycza_ce progu perkolacji modelu oraz korelacji kretosci z powierzchnia. 
charakterystyczna. i porowatoscia. efektywna, ukladu. 

5.1. Wprowadzenie 

Jak pokazalismy w poprzednim rozdziale, w przypadku procesow transportu przez 
osrodek porowaty, bardzo istotne jest, aby oprocz zjawisk dyfuzyjnych brae pod uwa- 
ge hydrodynamiczny przepb/w cieczy. Transport ten zachodzi w mikroskali przez 
system kanalikow, ktorych struktura i skomplikowana siec polqczen moze wplywac 
w sposob znaczqcy na jego wlasciwosci. Zwykle do opisu osrodkow porowatych uzy- 
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wa sie dwoch parametrow: porowatosci <p i przepuszczalnosci k Parametry te nie 
mowia_ jednak nic konkretnego o detalach mikrostruktury ukladu, ani o tym, jak za- 
chodzi transport. W tym kontekscie atrakcyjnym pomyslem wydaje sie zdefiniowa- 
nie dodatkowego, bezwymiarowego parametru, ktory opisywalby srednie wydluze- 
nie drogi czasteczek transportowanych przez osrodek porowaty, nazwanego w lite- 
raturze kretoscia_ {ang. tortuosity). Wtej czesci zajme sie badaniami numerycznymi 
nad kretosciq oraz zaleznoscia. kretosci od porowatosci w wybranym modelu osrod- 
ka porowatego. 

5.2. Model 

Model osrodka porowatego, dla ktorego liczyc bede kretosc, zdefiniowany jest na sie- 
ci kwadratowej o rozmiarze Lx Liu (starych sieci). Na sieci roztozone zostary losowo 
identyczne przeszkody w ksztalcie kwadratow o rozmiarach a* alu, gdzie 1 ^ a ^ L. 
Przeszkody te moga. dowolnie pokrywac sie nawzajem i sa_ calkowicie nieprzepusz- 
czalne dla cieczy Pozostab/obszar jest dla cieczydostepnyitu zachodzijej transport. 
Na rysunkuEUprzedstawione zostab/ uktady dla kilku wybranych porowatosci, przy 
czym kolorem szarym zaznaczono przeszkody reprezentuja.ce czesc nieprzepusz- 
czalnq osrodka porowatego. W celu minimalizacji wpb/wu efektow brzegowych na 
obliczenia, zalozylismy periodyczne warunki brzegowe w obu kierunkach (67l l68l . 
W celu wymuszenia przepb/wu cieczy przez mikrostrukture, na uklad nalozona zo- 
stala stala sila zewnetrzna (grawitacja) rownolegla do jednej z osi ukladu. 

5.2.1. Prog perkolacji 

Ze wzgledu na obliczenia kretosci w funkcji porowatosci ukladu, istotnym parame- 
trem jest upakowanie krytyczne, czyli prog perkolacji <p c . Okazuje sie, ze wielkosc 
przeszkod uzytych do budowy modelu ma znaczqcy wpb/w na polozenie progu per- 
kolacji. W naszych badaniach, podobnie jak w (67j[68]> we wszystkich symulacjach 
przyjeto a = 10. Przeprowadzone zostab/ obliczenia numeryczne progu perkolacji 
przypomocy algorytmunumerowaniaklastrowHoshena-Kopelmana 1 10 II orazme- 
tody Kirkpatricka wyznaczania punktu krytycznego 11021 dla a = 10 i otrzymalismy 
prog perkolacji (p c » 0.367. Wartosc ta lezy pomiedzy wartosciq <p ( P ~ 0.33(3) 11031 
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c) (p = 0.7 d) (p = 0.8 

Rysunek5.1. Przyklady konfiguracji na sieci o rozmiarach 200 x 200 {lu. x lu) (jednostek sie- 
ci) utworzonych przez losowe rozmieszczenie pokrywajqcych sie kwadratow 
o rozmiarach 10 x 10 [lu x lu) dla roznych porowatosci <p. 

(obliczona_ dla perkolacji ciqglej, a — ► oo), a wartoscia. (p^ ~ 0.5927 11041 (standar- 
dowa perkolacja wezlow, a = 1). Dodatkowo, na rysunku l5~2l przedstawiona zostala 
prosta analiza potozenia progu perkolacji z uzyciem metody skalowania rozmiaru. 
Na rysunku wyraznie widac, ze <p^ < <Pc < 0c? ■ Warto zauwazyc, ze obliczona przez 
nas wartosc <p c jest rozna od wartosci podawanych w (671 [68], gdzie autorzy blednie 
powiazali <p c dla a = 10 z <p^ ] dla a —* oo i uzywali (p c = (f>*p ~ 0.33(3). Dodatkowo na 
rysunku [531 przedstawiam zaleznosc progu perkolacji (p c od odwrotnosci dlugosci 
boku pojedynczej przeszkody a z uzyciem osi pollogarytmicznych. Wartosci nume- 
ryczne wyznaczone zostary algorytmem Hoshena-Kopelmana dla przeszkod wielko- 
sci 2 k , gdzie k = 0,1, 2... 6. Widac wyraznie, jak sukcesywne zwiekszanie rozmiaru 
przeszkody prowadzi do przesuniecia progu perkolacji od <pf ] dla klasycznej perko- 
lacji wezlow (ang. site percolation) do wartosci (p^ ] dla perkolacji cia_glej. 
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Rysunek 5.2. Jakosciowa analiza progu perkolacji dla a = 10 i trzech roznych rozmiarow sie- 
ci. Na osi y zaznaczone zostalo prawdopodobieristwo, ze uklad perkoluje. Na 
osi x odlozono porowatosc ukladu. Widac wyraznie, ze skalowanie wskazuje 
na prog perkolacji w okolicach (f> « 0.367 (punkt przeciecia krzywych). War- 
tosc (p^p ~ 0.33(3) obliczona zostala dla perkolacji ciaglej w 11031 , a wartosc 
(p^ ~ 0.5927 wyznaczono dla standardowej perkolacji wezlow 11041 . 
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Rysunek 5.3. Zaleznosc progu perkolacji cf> c od odwrotnosci dlugosci boku przeszkody 1/a. 

Symbole (o) reprezentujq numerycznq wartosc <p c w zaleznosci od rozmiaru 
przeszkody a = 2 lu (statych sieci), gdzie k — 0, 1,2. . . 6. Liniami przerywanymi 
zaznaczone zostaiy progi perkolacji ci^glej (p^ 1 oraz standardowej perkolacji 
sieciowej (f>f ] . 

5.3. Algorytmy numeryczne 

W celu wyznaczenia kretosci w badanym modelu musimy uzyc szeregu metod nu- 
merycznych. Poczawszy od generowania mikrostruktury, przez rozwia_zanie rownan 
przeptywu, wyznaczenie linii pra_du, wtasciwe obliczenia kretosci, az po weryflkacje 
czasu relaksacji w ukladzie, ekstrapolacje czesci rezultatow do stanow stacjonarnych 
oraz analize btedow numerycznych. 
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5.3.1. Konstrukcja osrodka porowatego 



Konstrukcja osrodka porowatego o okreslonej porowatosci jest realizowana przy po- 
mocy prostego algorytmu losowego osadzania swobodnie pokrywajqcych sie prze- 
szkod na badany obszar (63EH1- Rozpoczynajqc z pustego obszaru [cp = 1), nieprze- 
puszczalne kwadraty sq. na niego nakladane tak dtugo, jak dtugo (p > (po, gdzie (po jest 
oczekiwanq porowatosci^. 



5.3.2. Przepfyw cieczy 

Na potrzeby rozwiazania problemu przepfywu cieczy przez uklad porowaty zaadop- 
towany zostal model gazu sieciowego Boltzmanna (LBM) |84| w przyblizeniu BGK 
operatora kolizji |88| (patrz rozdziatO . Metodatapotwierdzitaswoj^uzytecznosc w 
symulacjach przeptywu cieczy w roznych warunkach (59111051 . 

Jednym z jej ograniczefi jest minimalna skala przestrzenna 4 lu, na ktorej jest 
mozliwe uzyskanie rozwiazania rownan przeprywu Naviera-Stokesa z zadowalajqcq 
dokladnosciq. Ograniczenie to ma dose duze znaczenie w przypadku porowatosci 
(p — - <pc, gdyz w tym przypadku znacznie rosnie liczba wqskich kanalow o szerokosci 
d < 4 lu. Z tego powodu na uklady skonstruowane wg opisanego powyzej algorytmu 
nalozona jest standardowa procedura zwiekszenia podzialu sieci na wiekszq ilosc 
komorek (rysunek l5.4l l. Procedura ta dzieli kazd^ z L 2 komorek sieci obliczeniowej 
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Rysunek 5.4. Procedura dodatkowego podzialu sieci osrodka porowatego. Stopien podzialu: 
a) fc re f = 2, b) k ie { = 4. Litery A, B i C uzyte zostaly do odroznienia odpowiada- 
jacych sobie elementow. 



na kwadraty Ar re f x A; re f, gdzie fc re f =1,2,... jest stopniem podzialu. Otrzymana siec ma 
zatem rozmiar A; re fL x A; re fL. 

Po inicjalizacji ukladu, petla obliczeniowa modelu jest kontynuowana f max kro- 
kow czasowych, gdzie t max jest wybrane wg kryterium, ktore omowimy w dalszej 
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czesci pracy (patrz podrozdzial l5.4.2ll . Dzigki uzyciu odpowiednich warunkow brze- 
gowych nakomorkach nieprzepuszczalnych {ang. midgrid) otrzymane rozwiazania 
majq dokladnosc drugiego rzgdu w czasie i przestrzeni |84 1 . Przyktadowe pola prg d- 
kosci uzyskane za pomocq omawianej metody dla kilku roznych wartosci porowato- 
sci przedstawione zostafy na rysunku l531 . 




c) = 0.80 d) = 0.95 



Rysunek 5.5. Predkosc cieczy u = |u| obliczona na sieci obliczeniowej o rozmiarach 600 x 600 
odpowiadajacej sieci fizycznej 200 lu x 200 lu (poziom podzialu /c re f = 3). Wiel- 
kosc przeszkod: 10 lu x 10 lu (t.j., 30 x 30 wezlow sieci). Porowatosc ukladu: 
a) (p = 0.45, b) (f> = 0.65, c) (f> = 0.8, d) (f> = 0.95. Warunki periodyczne zasto- 
sowane zostaly w obu kierunkach. Szare kwadraty reprezentujq obszar niedo- 
stepny dla cieczy, pozostaly obszar jest dozwolony dla przeplywu. Zewnetrzna 
sila grawitacyjna dziala w kierunku pionowym. 
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5.3.3. Linie pradu 

Kluczowq procedure calej operacji liczenia kretosci przeplywu jest dokladne wyzna- 
czenie linii prqdu (95). Aby je wyznaczyc, potrzebna jest znajomosc pola predkosci 
u(r) w calej objetosci osrodka. W naszym przypadku pole u(r) zostato uzyskane za 
pomocq interpolacji dwuliniowej {ang. bilinear interpolation) [ 1061 z predkosci we- 
zlowych u otrzymanych metoda. LBM. Na tak otrzymanym polu predkosci rozwiazy- 
wane jest rownanie ruchu bezmasowych czqstek (znacznikow): 

dr 

■TT=u(r), (5.1) 
at 

ktorych tory ruchu odpowiadajq, w przypadku stanu stacjonarnego, liniom prqdu 
przeph/waja_cej cieczy 1551 1 1071 . Ze wzgledu na skomplikowanq postac brzegow oraz 
znaczne roznice w lokalnych wartosciach interpolowanego pola predkosci u(r), do 
rozwia_zania rownari typu 115.111 zostala wykorzystana metoda Runge'go-Kutty czwar- 
tego rzedu ze zmiennym krokiem czasowym 11061 . 



5.3.4. Kretosc 

Zgodnie z rownaniem 11.21 . kretosc przeplywu zdefiniowana jako stosunek sredniej 
dlugosci linii prqdu przechodza_cych przez dany przekroj wjednostce czasu do sze- 
rokosci probki 1 1 1, co prowadzi do: 

I Uy{x)X{x)dx 

T =l—r ( 5 - 2 ) 

/ u y {x)dx 

J A 

gdzie A jest dowolnym przekrojem prostopadlym do osi y, a L oznacza dlugosc ukla- 
du. W rownaniu tym xe A, X(x) jest dlugosciq linii prqdu przecinajqcej A w punkcie 
x, a u y {x) - skladowa, wektora predkosci u w punkcie x prostopadlq do przekroju A. 
W literaturze calkowania w rownaniu 1 15.21 1 definiujqcym kretosc wykonywane byty 
albo za pomoca, metody Monte Carlo (MC) (67l l68l fT08l . albo poprzez kwadratury 
11091 . W metodzie MC dlugosci linii pra_du przechodza_cych przez losowo wybrane 
punkty na wybranym przekroju sa_ usredniane z uzyciem odpowiednich wag. Z kolei 
metoda kwadratur jest realizowana przez przyblizenie rownania 15.21 : 

£ Uy(Xj)A{Xj)AXj 

1 ; 

T*-^— , (5.3) 

L ^Uy{Xj)AXj 

j 
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gdzie przedziafy Ax 7 - = xj+i - xj sa_ dyskretnymi podziaiami przekroju. W ogolno- 
sci oba powyzsze podejscia do obliczenia calek w rownaniu 115.21 1 sa_ rownowazne, 
okazuje sie jednak, ze mogq bye one w rozny sposob implementowane. Na przyklad 
w (67111081 calkowanie Monte Carlo zostalo przeprowadzone na liniach prqdu prze- 
chodzqeych przez losowo wybrane wezly lezqce w obszarze porow osrodka, w (68] 
linie prqdu przecinaly dokladnie wszystkie komorki w obszarze porow, tymczasem 
w (6T1I109I za punkty startowe wybrane zostaly wezly brzegowe na wlocie cieczy do 
ukladu. Wszystkie wymienione podejscia la_czy jednorodny rozklad punktow starto- 
wych wybieranych losowo bqdz regularnie w obszarze calego osrodka. Sposob wy- 
brania komorek startowych dla linii prqdu moze miec ogromne znaczenie dla obli- 
czanych wartosci kretosci, szczegolnie ze wzgledu na selektywny charakter przeply- 
wu w obszarze niskich porowatosci, gdyz w tym przypadku ciecz wybiera kilka glow- 
nych kanatow przez ktore przeplywa niemal cala jej objetosc (rysunek [575l . Z tego 
powodu sumy w rownaniu 115.31 1 mogq zawierac wiele wyrazow o znikomym wprywie 
na wartosc kretosci. Aby ominqc ten problem, na rownanie 115.311 nalozylismy waru- 
nek stalosci strumienia pomiedzy dwiema sqsiadujqcymi liniami prqdu, co mozna 
wyrazic wzorem: 

u y {Xj)Axj = const. (5.4) 
Wartosci punktow Xj na przekroju wyznaczone zostaly ze wzoru rekurencyjnego: 

rxj l rL 

/ u y {x)dx = — / u y {x)dx, j = l,...,N, (5.5) 
Jx hl N Jo 

gdzie xo = 0. Zastosowanie warunku stalego strumienia pomiedzy sasiednimi liniami 
prqdu powoduje, ze wyrazenie 1 15.31 1 upraszcza sie do: 

11 « 

T "lNf l HXj) ' (5 - 6) 

gdzie AT jest liczba, linii prqdu wzieta_ do usredniania. We wzorze tym wszystkie wyra- 
zy sumy sa_ tego samego rzedu wielkosci. 

Ostatecznie, w celu obliczenia kretosci przeplywu przez osrodek porowaty o zna- 
nym kierunku makroskopowym przeplywu, wybierany jest przekroj prostopadly do 
tego kierunku. Nastepnie, z rownania 115.51 1 znajdowany jest rozklad punktow xj na 
przekroju, tak aby spelniony byl warunek 15 .41 . Linie pradu sa. wyznaczone przez 
sledzenie toru ruchu cza_stek (rownanie 15.111 o masie m = startujqeych w dwoch 
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przeciwlegfych kierunkach ze znalezionych punktow Xj na przekroju. Ruch czqstek 
sledzony jest tak dtugo, az dotrq one do brzegow ukladu x = lub x = L. Z tak wy- 
znaczonych linii prqdu obliczana jest srednia wyrazona rownaniem D5.6I I. Przyktad 
linii prqdu wyznaczonych wg omawianej procedury zostat przedstawiony na rysun- 
kuES 




c) = 0.80 d) = 0.95 



Rysunek 5.6. Linie pradu wyznaczone z warunku stalego strumienia (rownanie l5.4l l pomig- 
dzy sasiadujacymi liniami pradu dla tych samych ukladow jak na rysunku lS31 
Pozioma linia zaznaczona dla ukiadu o (p = 0.65 reprezentuje przekroj y = L/2, 
na ktorym znalezione zostaty punkty startowe xj dla linii pradu. Etykiety wska- 
zuja; 1) punkt koricowy niekompletnej linii pradu, 2) punkt startowy niekom- 
pletnej linii pradu, 3) „martwy" obszar osrodka, w ktorym nie nasUppuje prze- 
plyw cieczy. 
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Na tym rysunku zaznaczona zostala rowniez linia prqdu, ktora nie znajduje drogi 
przez caty obszar osrodka porowatego od brzegu x = do x = L. Moze sie to zdarzyc 
w przypadku linii pr^du, ktore przecinajq obszary o bardzo niskich predkosciach lo- 
kalnych. Przykladowa linia pra_du z rysunku [5T6l star tuje z punktu (2) o bardzo wyso- 
kiej wartosci predkosci lokalnej i probuje przejsc przez obszar (1), gdzie predkosci 
sq bardzo niskie. Linie, ktore nie znajduja_ drogi przez caty osrodek, nie sq brane pod 
uwage w sumach w rownaniu ( 15.31 1 . Dzieki zastosowaniu procedury uwzgledniajqcej 
warunek stalego strumienia pomiedzy sqsiadujqcymi liniami prqdu, „martwe" ob- 
szary osrodka porowatego s^ praktycznie calkowicie pomijane przy wyborze punk- 
tow startowych xj . 

5.3.5. Ekstrapolacja do stanow stacjonarnych 

Obliczenia metodq LBM prowadzone dla warunkow i wg procedury omowionej w po- 
przednich podrozdzialach prowadzq do uzyskania stanu stacjonarnego przeptywu. 
Dopiero stan asymptotyczny (stacjonarny) pozwala w efektywny sposob wyznaczyc 
linie prqdu potrzebne do obliczen kretosci. W celu wyznaczenia minimalnej ilosci 
krokow symulacji potrzebnych do otrzymania stanu stacjonarnego, badalem chwi- 
lowe wartosci T w kolejnych krokach czasowych frysunek 15.711 . Kazdy z punktow 

1.5 
1.4 
1.3 
1.2 
1.1 
1 

2500 5000 7500 10000 12500 15000 

t 

Rysunek 5.7. Zaleznosc T od kroku czasowego w procesie relaksacji do stanu stacjonarnego 
dla kilku roznych porowatosci ukladu (p. Kazdy z punktow przedstawionych 
na wykresie stanowi wynik usredniania po 25 ukladach. Liniami zaznaczone 
zostaly dopasowane relacje wykladnicze postaci 15.71 . 

pomiarowych przedstawionych na rysunku ISTTl jest kretosciq usrednionq po N = 25 



4>=0.55 — 

4>=0.65 >--■*> 
$=0.75 

(j>=0.85 ----- 
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roznych konfiguracjach o okreslonej porowatosci (p. Okazalo sie, ze po t krokach, 
zaleznosc T(t) moze bye przyblizona relacjq wyktadniczaj 

T{t) a r s -cexp(-z7f r ei), t> t , (5.7) 

gdzie T s , c oraz t ie \ sq pewnymi parametrami i rozniq sie wartosciami dla roznych 
porowatosci 0. Z rysunku widac wyraznie, ze czas relaksacji uktadu rosnie wraz ze 
zblizaniem sie porowatosci (p do 1, co przedyskutuje dokladniej w dalszej czesci pra- 

5.4. Wyniki 

5.4.1. Wartosci lokalne na przekroju 

Wartosci wielkosci potrzebnych do wyznaczenia calek w rownaniu 15.2H zostafy przed- 
stawione na rysunkach l5.8l - [5TT0l Na przekrojach poziomych obu ukladow obliczono: 

1. znormalizowanq predkosc gdzie u^^ix) jest najwiekszq pred- 
kosciq wyznaczona, wzdluz wszystkich linii pradu w uktadzie, 

2. kretosc pojedynczej linii prqdu t(x) = X[x)/L, 

3 . ilo czyn t (x) u y (x) / w™ 3 * (x) , 

4. stosunek minimalnej do maksymalnej predkosci wzdluz linii pra_du. 

Wszystkie wielkosci zostary zebrane z pol predkosci dla ukladow z rysunkow l5.5l oraz 
15.61 Zaleznosci przedstawione na wykresach Il5.8l45.10l a-c) sq silnie zalezne od po- 
zycji, na ktorej byl wyznaczony przekroj. W naszym przypadku przekroje byry wziete 
dla y = L/2. Jak mozna bylo sie spodziewac, profile predkosci u y {x)l u^^ix) s^ ciq- 
gle i kawalkami rozniczkowalne, czego nie mozna powiedziec o przebiegach kretosci 
t(x). Co oczywiste, iloczyn T(x)u y (x)/u™ ax (x) rowniez nie jest ciqgry i jest to jeden 
z powodow trudnosci, jakie pojawiajq sie przy calkowaniu numerycznym rownania 
D5.2I I. Kazda nieciqglosc funkeji t(x) przedstawionej na omawianych wykresach w 
czesci b) odpowiada laczeniu sie lub rozwidlaniu dwoch sasiadujacych linii pradu. 
Odpowiada to sytuacji, kiedy dwie strugi cieczy napotykajqee przeszkode oph/wajajq 
z dwoch roznych stron. Poprzez liczbe punktow nieciqglosci funkeji t(x) oszacowac 
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Rysunek 5.8. Lokalne wartosci obliczane na przekroju A uzywane do obliczeri kretosci: (a) 
skladowa predkosci u y [x) znormalizowana do wartosci maksymalnej u™ 3 * w 
5 x 10" 5 lu tu~ l ; (b) lokalna krgtosc t{x) = \{x)IL; (c) iloczyn r(x)u j ,jc/u™ ax 
i (d) stosunek prgdkosci minimalnej do maksymalnej (obie wartosci znalezio- 
ne wzdluz calej linii pradu). Wartosci z tego rysunku zostaly zebrane z prze- 
kroju dla konfiguracji o porowatosci (p = 0.65 z rysunku lSTSl 



mozna liczbe tego typu „wysp" w ukladzie. Porownujac profile dla roznych porowa- 
tosci od cf> = 0.45 do (p = 0.9 mozna zauwazyc, ze liczba punktow niecia_glosci maleje 
dla (f>^ (p c oraz dla 0—1. Problem znajdowania miejsc niecia.gk>sci funkcji t{x) na 
przekroju nie jest dobrze okreslony numerycznie. Wybor miejsc startowych dla linii 
prqdu blisko punktow niecia_glosci powoduje, ze wartosci kretosci moga. bye obar- 
czone duzymi bledami, gdyz niewielka niedokladnosc w wyznaczeniu xj moze skut- 
kowac duzq niedokladnoscia. (skokiem) wartosci X(Xj). Okazuje sie rowniez, ze sto- 
sunek predkosci maksymalnej do minimalnej rosnie znaczqco w punktach bliskich 
punktom niecia_gk>sci. Jednym z czynnikow wplywaja_cych na dokladnosc zwiaza- 
nq z nieciqglosciami w t(x) moze bye liczba linii pra_du wzietych do obliczen. Aby 
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x 

Rysunek 5.9. Lokalne wartosci obliczane na przekroju y = L/2 uzywane do obliczeri krgtosci 
dla konfiguracji o porowatosci cf> = 0.80 z rysunku [5?6l Pozostale parametry jak 
na rysunku [5^8l 

sprawdzic wpbyw tego parametru na otrzymywane wartosci, dla kilku ukladow o roz- 
miarach 200 x 200 lu (stafych sieci) z podzialem k ie f = 3 o roznych (p wyznaczona 
zostala wartosc T obliczona zN = 2 k linii pra_du, gdzie k = 1, 2, . . . , 10. Wyniki przed- 
stawione zostary na rysunku [5.11l Widac wyraznie, ze w okolicach N « 100 krgtosc 
przestaje sie- zmieniac wraz z zage-szczaniem linii prqdu. Mozna stqd wnioskowac, ze 
liczba N « L linii prqdu jest wystarczajqca do obliczen kretosci, a jej zwiekszenie nie 
zmieniloby znacz^co otrzymywanych wartosci T. 



5.4.2. Czas relaksacji t m \ 

Jak pokazalismy w podrozdziale l5.3.5l f rysunek \577\ wartosc T jest zalezna od kroku 
czasowego t symulacji. Charakterystyczny czas relaksacji tego zjawiska t Te \ jest roz- 
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Rysunek 5.10. Lokalne wartosci obliczane na przekroju y = L/2 uzywane do obliczeri kreto- 
sci dla konfiguracji o porowatosci (p = 0.95 z rysunku lSTSl Pozostale parametry 
jak na rysunkuEEJ 



1.8 
1.7 
1.6 
1.5 
I- 1.4 
1.3 
1.2 
1.1 
1 





• 










A 










A A 


A- - 








___A--- 


- - - A 








- - -A - - . 


- - - -A - " ■ 














fl\ Q 








— e — 


— e — 


— G 
















































— v- - 


— v~ 


— v— 


— V 






▼ T - - ■ 










- -T- - 


- - ▼ 



10 100 
N 



1000 



=0.45 — e- 

=0.50 • 

=0.55 •• 

=0.60 - 
=0.70 

=0.90 - 



Rysunek 5.1 1. Kretosc w funkcji liczby linii prqdu dla ukladu 200 x 200 lu i podzialu fc re f = 3. 
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ny dla roznych porowatosci oraz roznych konfiguracji porow w osrodku. Zaleznosc 
f re i(0) przedstawiona zostata na rysunku l5.121 gdzie kazdy punkt reprezentuje sred- 
niq po co najmniej 25 ukladach. Jak widac z wykresu, czas relaksacji ma minimum 
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Rysunek5.12. Czas relaksacji £ re i w funkcji porowatosci </> dla L = 200 Zu, /c re f = 3 i f max = 
30000 tu. 



dla « 0.6 i od tej wartosci rosnie w obu kierunkach. Wzrost w kierunku wysokich 
porowatosci {(p — - 1) jest jasny i zgadza sie z innymi doniesieniami o takim zacho- 
waniu modelu LBM w zakresie niskich liczb Macha 1 1 10 1 . Wzrost f re i w kierunku wy- 
sokich porowatosci {(p — ► (p c ~ 0.367) jest prawdopodobnie zwiazany ze wzrostem 
liczby i dlugosci martwych pol w osrodku porowatym, ktore sq wypelnione cieczq, 



ale praktycznie nie dajq zadnego wkladu do przepfywu 11111 . Czysto dyfuzyjny cha- 
rakter transportu w tych porach powoduje znaczqcy wzrost f re i. 

Jak pokazalismy w podrozdziale 15.3.51 w zaleznosci T{t) wystepuje czynnik to- 
Okazuje sie, ze czynnik ten spelnia warunek 500 < to < 1000 w cah/m zakresie (p. 
Dlatego, korzystajqc z danych z rysunku [5.121 do obliczen kretosci zalozylismy mak- 
symalnq liczbe f max = 1.5 x lO 4 krokowdla0 < 0.8 i3x 10 4 dla0 > 0.8. Okazalo siejed- 
nak, ze dla kilku uktadow z porowatosci^ (p = 0.45 to bylo znacznie wieksze od 10 3 tu. 
Jesli przyjrzec sie blizej strukturze tych ukladow, zauwazyc mozna, ze sq one zbudo- 
wane w odmienny sposob. Posiadajq one dwa lub wiecej glownych kanalow, ktorymi 
zachodzi transport wzdluz sih/ wymuszajqcej przeph/w i sq one polqczone ze sobq 
kanalem mniej wiecej prostopadtym do kierunku przeph/wu glownego. W ukladzie 
takim, w procesie relaksacji do stanu stacjonarnego dojsc moze do „przelqczenia" 



51 



przepfywu z rownolegfych kanalow do kanalu la_cza_cego, co widoczne jest na wykre- 
sie T(t) jako nagty uskok wartosci T. Przyktad relaksacji kretosci dla przykladowego 
ukladu o porowatosci przedstawiony zostat na wykresie l5.131 Ze wstawionych do 




1 1 * 1 ' ^ 

5000 10000 15000 20000 

t 

Rysunek5.13. Przyktad zbieznosci wartosci T do stanu asymptotycznego. Na rysunkach 
wstawionych do wykresu widac konfiguracje linii prqdu odpowiadajace cza- 
som t = 2500,4500 i 17000. Sila zewnetrzna skierowana jest w kierunku pio- 
nowym. 

wykresu wizualizacji linii prqdu dla trzech roznych chwil czasu widac wyraznie, ze 
wybrany do przeph/wu kanat pionowy, po okolo 4000 krokow przestaje bye jedynym 
dozwolonym dla przeph/wu, a wokolicach t = 15000 przestaje dominowac. Co cieka- 
we, przez wybranie poczadcowo najkrotszej sciezki, uklad wydawal sie bye stabilny 
ze wzgledu na kretosc az do okolo t = 4000, a nastepnie z powodu znalezienia do- 
datkowej drogi dla przeph/wu, wartosc T zaczela gwaltownie rosnqc, by nastepnie 
kontynuowac wzrost z relaksacja. wykladnicza,. Jak widac z powyzszej analizy, w te- 
go typu ukladach trudno jest zdefiniowac scisle kryterium zbieznosci wartosci T do 
stanu asymptotycznego. Dlatego nie ma pewnosci, ze we wszystkich badanych ukla- 
dach w okolicy (p = 0.45 osiajmiety zostat stan stacjonarny. Wartosci T(0 = 0.45) mo- 
ga_ bye wiec obarczone niewielkim, systematycznym btedem, ktorego wielkosc osza- 
cowac mozna na mniej niz 1% (na podstawie liczby ukladow, ktore zdradzaja. wyzej 
wymienione cechy). Minimalna wartosc czasu relaksacji w okolicach (p « 0.6 moze 
bye zwia_zana z przeprowadzonq w poprzednim podrozdziale analiza_ ilosci punktow 
niecia_glosci w zaleznosci od (p. W okolicach (p « 0.6 niecia_gk>sci jest bardzo duzo 
(najwiecej) co odpowiadac moze selekcji kanalow o niskiej kretosci, ktorych ilosc jest 
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dose duza (kazdy punkt nieciqglosci to taki wlasnie kanal) . Mozna sie spodziewac, ze 
relaksacja wkanatach o takiej charakterystyce bedzie najszybsza. 



5.4.3. Znaczenie wielkosci ukfadu 

Jednq z wazniejszych cech kretosci ukladu, jakq udalo mi sie zauwazyc, a jaka potrak- 
towana zostala pobieznie w poprzednich pracach innych autorow, jest jej zaleznosc 
od wielkosci sieci L. W celu ilosciowej analizy tej zaleznosci, srednia kretosc T obli- 
czona zostala dla trzech porowatosci: (p - 0.5, 0.7 i 0.9 oraz czterech wielkosci ukladu 
L = 50, 100, 200 i 300. Wyniki przedstawione zostaly na rysunku [5.141 Kazdy z punk- 
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Rysunek5.14. Zaleznosc T{L) dla c6 = 0.5, 0.7 i 0.9. Punkty reprezentujq wartosci srednie, 
usrednione po N ^ 24 ukladach. Krzywe dopasowaniami do wzoru 15.8L a 
slupki bledow reprezentujq blad standardowy wartosci sredniej. 
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tow wykresu reprezentuje sredniq wartosc T usrednionq. po N ^ 24 ukladach (./V do- 
brane zostalo tak, by bledy statystyczne widoczny na wykresie byby tego samego rze- 
du). Krzywe widoczne na wykresie sq dopasowanymi zaleznosciami wykladniczymi 
wyrazonymi wzorem: 

T[L) = T 00 -bexp[-cL), (5.8) 

gdzie jest wartosciq asymptotycznq kretosci (dla ukladu nieskonczonego) i razem 
z b oraz c tworzy grupe parametrow swobodnych dopasowania. Mozna zauwazyc, 
ze w cafym zakresie porowatosci T jest rosnaca^ funkejq rozmiaru ukladu. Co wiecej, 
wielkosc charakterystyczna L* powyzej ktorej T nie zmienia sie znaczqoo wraz ze 
wzrostem L wynosi L* « 200. Ponadto, rysunek r5.14l nie pozostawia wqtpliwosci, ze 
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zignorowanie efektow zwiazanych z rozmiarem sieci prowadzi do niedoszacowania 
wartosci T. Wielkosc niedoszacowania jest zalezna od porowatosci (p i spodziewac 
sie mozna, ze wieksze roznice wysta_pic moga_ dla nizszych porowatosci. 

Kolejnym krokiem bylo sprawdzenie, jak na otrzymywane wartosci wptywa sto- 
pien podziaiu sieci fc re f. Na rysunku [5.151 przedstawione zostary wartosci T w funk- 
cji wielkosci ukladu dla (p = 0.5 oraz czterech stopni podziaiu sieci fc re f = 1,2,3 i 4. 
Okazato sie, ze mimo znanego z literatury faktu, ze model LBM potrzebuje przynaj- 
mniej cztery jednostki sieci do poprawnego opisania hydrodynamiki (84) (co w na- 
szych warunkach odpowiada podzialowi /c re f = 3), wprzypadku obliczen kretosci sto- 
pien podziaiu nie ma tak duzego znaczenia. Z wykresu ( 15. 151 1 widac wyraznie, ze juz 
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Rysunek 5.15. Zaleznosc T[L] dla cp = 0.5 i czterech stopni podziaiu fc re f = 1,2,3 i 4. Punkty 
reprezentujq wartosci srednie, usrednione po n = 9 ukladach. Slupki bledow 
reprezentujq blqd standardowy wartosci sredniej. 

fc re f = 1 wydaje sif bye wystarczaja.ee do tego typu obliczen. Ten niespodziewany wy- 
nik (czesc ukladu na przewezeniach posiada przeciez lokalne rozwiqzania niezgodne 
z rownaniami Naviera-Stokesa) mozna wytlumaczyc zauwazajqc, ze dla porowato- 
sci duzo wyzszych od progu perkolacji glowny przeplyw cieczy przez uklad zacho- 
dzi szerokimi kanalami o dose duzej srednicy w porownaniu z komorkq elementar- 
nq sieci. Dlatego wplyw niedokladnosci w wyznaczeniu pola predkosci w okolicach 
przewezeri nie ma istotniejszego znaczenia. Spodziewac sie jednak mozna, ze wraz 
ze zblizaniem sie z porowatosci^ do progu perkolacji [<p —> (p c ), wystepowanie szero- 
kich kanalow dominujqeych przeplyw bedzie rzadsze, a co za tym idzie - znaczenie 
stopnia podziaiu moze bye wieksze. 
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5.4.4. Zaleznosc kretosci od porowatosci 

Dzieki znajomosci minimalnych wymagan odnosnie rozmiaru sieci, stopnia podzia- 
lu oraz dlugosci czasu relaksacji do stanu stacjonarnego, jestesmy w stanie wyzna- 
czyc zaleznosc kretosci T od porowatosci w szerokim zakresie (p. Dla porowato- 
sci <p = 0.45, 0.5,..., 0.95 wybrane zostafy: L = L* oraz stopien podzialu k ie f = 3. Dla 
kazdej porowatosci, kretosc T wyznaczona zostala dla M roznych konfiguracji prze- 
szkod, gdzie M zmienialo sie od 25 dla (p = 0.95 do 100 dla (p = 0.45. Srednie war- 
tosci T przedstawione zostary na wykresie 15.1 61 wraz z relacja. wyrazona. rownaniem 
H1.8el l dla tego samego zagadnienia zaproponowana. przez Koponena i innych w (68) . 
Przyczyny roznic pomiedzy wynikami prezentowanymi tutaj a praca. (68| mozna tlu- 
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Rysunek5.16. Zaleznosc kretosci T od porowatosci (p. Dane otrzymane zostary wg omo- 
wionej procedury z uzyciem modelu LBM i rownan l !5.6t oraz H5.7H (krzyze 
ze slupkami bledow). Dodatkowo zaznaczone zostary: 1) krzywa teoretyczna 
zaproponowana dla takiego samego modelu przez Koponena i in. |68| (linia 
przerywana) oraz najlepsze dopasowanie do rownania l ll.8bt (linia ci^gla). 

maczyc nieuwzglednieniem przez autorow tamtej pracy efektow wielkosci sieci, wy- 
borem L< L* oraz brakiem dokladnej analizy czasow relaksacji kretosci. 

Dane z rysunku [5 . 1 61 dopasowane zostaby do czterech jednoparametrowych wzo- 
row empirycznych z jednym wolnym parametrem zaproponowanych w innych opra- 
cowaniach (rozdzial ll.31 str. W przypadku naszych obliczen kretosci wyznaczo- 
nej dla modelu pokrywaja_cych sie prostokqtow przy pomocy procedury omowio- 
nej powyzej, najlepsze dopasowanie dal wzor empiryczny Jl.8bl l z parametrem p = 
0.77 + 0.03. Najlepsze dopasowanie zostalo narysowane liniq cia_gta_ na rysunku l5.161 
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Okazuje sie, ze dla (p « 1 dane odstepujq nieznacznie od wzoru [Obi Niezbyt do- 
bre dopasowanie danych w tym zakresie porowatosci moze bye m.in. rezultatem 
nieuwzglednienia anizotropii badanych uktadow spowodowanej periodycznymi wa- 
runkami brzegowymi. 



5.4.5. Korelacja kr^tosci z powierzehnia charakterystyczna 

Fakt, ze najlepszym dopasowaniem do danych numerycznych zaleznosci T{<p) oka- 
zalo sie prawo logarytmiczne, ma dose ciekawe konsekwencje dla korelacji tej wiel- 
kosci z innymi parametrami makroskopowymi opisujqcymi strukture osrodka. Wia- 
domo bowiem, ze powierzehnia charakterystyczna S osrodka zlozonego z losowo 
rozmieszczonych, pokrywajqeych sie kwadratowych przeszkod o polu powierzeh- 
ni Vo i obwodzie A , charakteryzuje sie logarytmicznq zaleznosci^ od porowatosci 
(1681): 

d 

(5.9) 



S = --<p\n<p, 

gdzie R jest tzw. promieniem hydraulicznym przeszkod, a d jest wymiarem prze- 
strzeni (tu d = 2). Promien hydrauliczny definiujemyjako R = d- Vol Aq, gdzie Vb jest 
objetosciq przeszkody, a A jej powierzehniq. Na rysunku l5.17l przedstawione zosta- 
lo porownanie powyzszej relacji z wartosciami numerycznymi. Punkty na wykresie 
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Rysunek5.17. Zaleznosc powierzehni charakterystycznej S od porowatosci dla kwadratow 
10 x 10 i k re { = 3. Linia. przerywana_ zaznaczony zostal wynik analityczny i5M 
z uwzglednieniem wartosci d = 2 oraz R = 15. 



reprezentujq srednie po N = 14 ukladach dla kazdej z przedstawionych wartosci (p. 
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Ogolna zgodnosc jest bardzo dobra, a niewielkie odstepstwa od wartosci analitycz- 
nej w okolicy progu perkolacji mogq bye zwiazane z efektami rozmiaru sieci oraz 
uzyciem fc re f = 3. 

Przy pomocy prostych rachunkow, mozna pokazac, ze z relacji l ll.8bl l oraz 1 15.91 1 
wynika, ze 11121 : 

S 

r-loc -. (5.10) 

(p 

Na rysunku f5. 18l przedstawiony zostat wykres zaleznosci T— 1 wfunkcji 4, gdzie stala 
proporcjonalnosci p = 0.77 wyznaczona zostata z procedury dopasowania do wzoru 
<1.8bl . Okazuje sie, ze proporcjonalnosc wyrazona wzorem H5.10I I dobrze opisuje za- 

1.14 i , , . . , . , 1 
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Rysunek5.18. Graficzna weryfikacja proporcjonalnosci wyrazonej wzorem f5 . 1 1 Linia pro- 
sta reprezentuje dopasowanie liniowe postaci y = a ■ x + b, przy czym najlep- 
sze dopasowanie uzyskane zostaio dla a = 0.16 i b = 1.0. 

leznosc miedzy T, S i (p w cafym badanym zakresie porowatosci i kre tosci. Oznacza 
to, ze w badanym modelu wzrost dlugosci linii brzegowej (powierzehni charaktery- 
stycznej) S prowadzi do wzrostu kretosci T, a efekt ten wzmacnia sie wraz ze zmniej- 
szaniem porowatosci [cp w mianowniku) . Na tym etapie nie mozna jednak stwierdzic, 
czy relacja H5.10I I jest w jakimkolwiek stopniu ogolna. Warto jednak zauwazyc, ze na- 
wet jesli relacja pomiedzy kretosciq a porowatosci^ w jakims ukladzie nie jest loga- 
rytmiczna, to zaleznosc 115.101 1 moze bye nadal spelniona. Jest to jeden z kierunkow 
dalszych badan nad korelacjami parametrow makroskopowych w osrodkach poro- 
watych. 
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5.4.6. Korelacja kr^tosci z porowatoscia efektywnq 

Analize korelacji porowatosci z kretosciq mozna tez przeprowadzic dla porowatosci 
efektywnej e ff> ktorq tu definiujemy jako stosunek ilosci wezlow sieci, przez ktore 
zachodzi efektywnie caly transport cieczy, do reszty wezlow (68) • Dzieki tej definicji 
wszelkie strefy „martwe", ktore z punktu widzenia transportu hydrodynamicznego 
pomijalne (ale wypelnione cieczy, dlatego wchodzqee do wartosci porowatosci kla- 
sycznej), nie sq wliczane do wartosci eff . Wartosc eff jest wiec z zalozenia mniej- 
sza niz klasyczna porowatosc zdefiniowana jako stosunek objetosci porow dostep- 
nej dla cieczy do objetosci zajmowanej przez caly material porowaty. W celu nume- 
rycznego wyznaczenia porowatosci efektywnej, zastosowatem procedure przedsta- 
wionq w |68] uwzgledniajqc wprowadzonq w niniejszej rozprawie zmodyfikowanq 
procedure wyznaczania linii prqdu. Za komorki sieci, ktore aktywnie biorq udzial 
w transporcie cieczy, brane byly te komorki, przez ktore przechodzity linie pr^du 
wyznaczone metodq calkowania znacznikow o zerowej masie unoszonych przez po- 
le predkosci. Wyznaczona w ten sposob zaleznosc porowatosci efektywnej dla kilku 
roznych porowatosci przedstawiona zostata na rysunku [5.19[ W tym miejscu warto 
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Rysunek5.19. Zaleznosc porowatosci efektywnej (/> e ff od porowatosci calkowitej <p. Czarne 
punkty to rezultat procedury numerycznej. Dodatkowo narysowane zostafy 
funkcje y = x oraz (/> e ff = 1 + ln((/>) - postulowana zaleznosc </> e ff ((/>). 

zauwazyc, jakie konsekwencje mialo uwzglednienie blednej wartosci progu perko- 
lacji (p c przez autorow pracy (68) (patrz podrozdzial l5.2.1L ktorzy tez zauwazyli, ze 
zaleznosc (p e &(4>) zachowuje sie jak zaleznosc logarytmiczna. Ze wzgledu na przyje- 
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cie blednego (p c , nie byli jednak w stanie wyjasnic, dlaczego proponowana zaleznosc 
0eff = 1 _ ln(0)/ln(0 c ) daje wartosci niefizyczne (tzn. <£ e ff > (p) dla wysokich poro- 
watosci. Jesli zauwazyc, ze prawdziwe (p c przyjmowato wartosc wieksza^ od przyje- 
tej przez autorow (68) - problem ten znika. Dla kwadratow o boku a = 10, wiel- 
kosc ln(0 c ) « -1.0 i powyzsze wyrazenie na e ff(0) przyjmuje uproszczonq postac: 
0eff = 1 + ln(0). Widac wyraznie frysunek !5.19ll . ze przy uwzglednieniu poprawnej 
wartosci progu perkolacji, zaleznosc ta (dla tego modelu oraz tych wielkosci prze- 
szkod) jest logarytmiczna. Ma to bardzo interesujqce konsekwencje, bo rowniez lo- 
garytmiczne zaleznosci powierzchni swobodnej S i kretosci T pozwalajq postulowac 
inne zwiqzki, z ktorych najciekawszym wydaje sie proporcjonalnosc miedzy kreto- 
scia. a porowatosciq efektywna; 

TcX(/)eff. (5.11) 

Porownanie wzoru 15.111 z rezultatami symulacji zostalo przedstawione na rysun- 
ku !5.201 Widac wyraznie, ze w carym zakresie porowatosci efektywnej kretosc jest do 
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Rysunek 5.20. Liniowa zaleznosc kretosci T od porowatosci efektywnej (/> e ff. 

niej wprost proporcjonalna. Wspolczynnik nachylenia prostej wynosi a = -0.68. Na 
tym etapie nie ma jednak podstaw do interpretacji tej wartosci. Podobnie jak w przy- 
padku powierzchni swobodnej, przedstawiona_ proporcjonalnosc nalezatoby zbadac 
w szerszym zakresie ukladow o roznej budowie, co pozwoliloby potwierdzic ba.dz 
odrzucic ja_jako prawo uniwersalne. 
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5.5. Omowienie wynikow i podsumowanie 



Przedstawiona w tym rozdziale dokladna analiza procedury numerycznej do obli- 
czania kretosci oraz wyniki, jakie uzyskalismy, pokazujq, z jak trudnym problemem 
mamy do czynienia. Skomplikowana struktura linii pra_du, ktorych uzywalismy do 
obliczen kretosci i duze roznice predkosci minimalnych i maksymalnych wzdluz po- 
jedynczej linii spowodowaly potrzebe uzycia wyrafinowanych metod caikowania ze 
zmiennym krokiem czasowym. Okazalo sie, ze nawet z uzyciem takich metod nu- 
merycznych, linie prqdu nie zawsze znajduja_ droge przez caly uklad, moga_ rozszcze- 
piac sie na charakterystyczne „wyspy", a nastepnie trafic w obszary o bardzo niskich 
predkosciach lokalnych, gdzie nawet zastosowanie zmiennego kroku czasowego nie 
umozliwia uzyskania zadowalajqcego rozwiazania. Oszacowanie tych bledow poka- 
zalo jednak, ze jest to efekt, ktory nie wph/wa znaczqco na obliczana_ kretosc. 

Kolejnym krokiem w strone usprawnienia procedury numerycznego wyznacza- 
nia kretosci ukladu bylo wprowadzenie rozktadu linii prqdu uwzgledniaja_cego waru- 
nek stalego strumienia cieczy pomiedzy sqsiednimi liniami prqdu. Dzieki temu linie 
praxlu w obszarach o wiekszej predkosci zageszczone, a srednia z linii prqdu jest 
juz zwyklq sredniq arytmetycznq. Nie chodzi tu tylko o uproszczenie wzoru na kre- 
tosc, ale przede wszystkim o zwiekszenie liczby linii prqdu w miejscach szczegolnie 
waznych, a takimi sa. wlasnie miejsca przez ktore nastepuje przepryw. Miejsca o lo- 
kalnie niskich predkosciach nie wnoszq wiele do transportu i nie sa_ brane pod uwage 
przy obliczeniach. Realizowana wczesniej w literaturze procedura regularnego roz- 
kiadu linii pra_du w pola_czeniu z obliczaniem sredniej wazonej z wagami rownymi 
predkosciom lokalnym nie jest w stanie dac dokladnosci, jakq oferuje algorytm za- 
proponowany w niniejszej rozprawie. 

W trakcie obliczen okazalo sie rowniez, ze metoda gazu sieciowego LBM jest bar- 
dzo dobrym narzedziem dla srednich porowatosci. Dla (p — ► (p c oraz 0^1 poja- 
wiaty sie rozne, niepozqdane efekty. Przede wszystkim, dla obu skrajnych wartosci 
(p znacznie rosnie czas relaksacji f re i kretosci. Oznacza to koniecznosc znacznego 
zwiekszenia liczby krokow czasowych potrzebnych do uzyskania stanu stacjonarne- 
go oraz - w niektorych przypadkach - ekstrapolacji tych wartosci za pomocq funk- 
cji wykladniczej. Zrodlo zwiekszenia czasu relaksacji jest jednak w obu przypadkach 
inne. Dla (p « 1 uklad zbudowany jest z pojedynczych przeszkod, ktore nie tworza_ 
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skomplikowanej sieci kanalow, a ich liczba jest tak mala, ze stan stacjonarny osia_- 
gany jest bardzo powoli (czego powodem jest dyfuzyjny charakter relaksacji pola ci- 
snienia). Jest to znana cecha modeli opartych o gaz sieciowy. W drugim skrajnym 
obszarze, bliskim progowi perkolacji, wzrost r re i mozna tlumaczyc istnieniem sie- 
ci dtugich, wzajemnie pol^czonych kanalow o malej srednicy, w ktorych relaksacja 
cisnienia nastepuje rowniez w sposob dyfuzyjny. Okazalo sie, ze konfiguracje o po- 
rowatosci w okolicach (p = 0.6 maja_ budowe optymalnq dla czasu relaksacji. Jest to 
obszar porowatosci w ktorym uklady maja_ charakter posredni - juz nie skladaja. sie 
jedynie z odizolowanych wysp, a jeszcze nie tworzq sieci dlugich, pojedynczych, od- 
izolowanych kanalow. 

Do analizowanej w tym rozdziale empirycznej relacji logarytmicznej pomiedzy 
kretoscia., a porowatoscia_ Jl .8bl> nalezy podejsc krytycznie. Dopasowanie, jakie zo- 
stalo uzyskane jest bardzo dobre, ale nie mozna zapominac, ze uzyskano je tylko 
dla modelu pokrywaja.cych sie przeszkod kwadratowych. Na tym etapie nie mozna 
wnioskowac, ze prawo to spelnione jest dla dowolnego ukladu, ani ze ma ono taka_ 
samq postac w przestrzeni trojwymiarowej. Dose obiecuja_co wygla_daja_ natomiast 
relacje, jakie udalo sie zaproponowac pomiedzy kretoscia. a charakterystyczna. po- 
wierzehniq mikroporow oraz efektywnq porowatoscia., gdzie uzyskana zostata bar- 
dzo dobra zgodnosc. Jednym z mozliwych dalszych kierunkow badan nad omawia- 
nymi zagadnieniami jest dalsze poszukiwanie zwiqzkow miedzy tymi parametrami 
makroskopowymi, zarowno dla roznych modeli, jak i w badaniach eksperymental- 
nych. 
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Rozdziaf 

6 



Podsumowanie 



Gtownym celem niniejszej rozprawy bylo zbadanie dwoch szczegolnych zagadnien 
zwiazanych z transportem plynow przez osrodki porowate: dyfuzji „anomalnej" wmo- 
delu Kiintza i Lavallee'go oraz analiza kretosci przeph/wu przez osrodek porowaty. 

W rozdzialach|2]i|4]wprowadzony zostat model Kiintza i Lavallee'go (KL) osrodka 
porowatego zlozonego z losowo rozmieszczonych rozpraszaczy. W niniejszej rozpra- 
wie udalo sie pokazac, ze obserwowana w tym modelu ewolucja profilu koncentracji, 
uprzednio interpretowana jako efekt dyfuzji anomalnej, moze bye wyjasniona jako 
zlozenie dwoch klasycznych zjawisk: dyfuzji Ficka oraz stabego dryfu hydrodyna- 
micznego 11131 . 

Naturalnq konsekwencjq tych badan bylo przejscie do badan nad zjawiskami hy- 
drodynamicznymi w osrodkach o bardziej skomplikowanej budowie. Problemem, 
ktory l^czy oba zagadnienia ze sobq, jest zbadanie korelacji pomiedzy makroskopo- 
wymi wielkosciami fizycznymi opisujqcymi transport cieczy przez osrodek porowa- 
ty. Po wprowadzeniu modelu gazu sieciowego Boltzmanna w rozdziale|3j omowio- 
na zostala nowa procedura numeryczna do wyliczenia kretosci (T) - parametru ma- 
kroskopowego opisujqeego stopiefi wydtuzenia drogi transportowanych przez osro- 
dek czqsteczek materii (rozdzialO. Dzieki dokladnej analizie bledow o roznym zro- 
dle pochodzenia, udalo sie przedstawic nowe wyniki dotyczqee zaleznosci pomie- 
dzy porowatoscia_ a kretosciq osrodka 11121 . Nowe, dokladniejsze wyniki dotyczqee 
zaleznosci kretosci od porowatosci pozwolify wysunac hipoteze dotyczqeq korelacji 
pomiedzy takimi parametrami jakporowatosc efektywna, powierzehnia charaktery- 
styczna i kretosc. 
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Kolejnym krokiem w kontynuacji przedstawionych w niniejszej rozprawie badan 
bedzie sprawdzenie uniwersalnosci przedstawionych korelacji. Rownie istotne wy- 
daje sie bye sprawdzenie, jak badane wielkosci zachowujq sie w uktadach (mode- 
lach) trojwymiarowych. Bardzo wazne, szczegolnie z punktu widzenia ewentualnych 
zastosowan, byioby rowniez opracowanie dokiadnych i efektywnych metod pomiaru 
eksperymentalnego wielkosci takich jak e ff, S i T. Mogh/by w nich znalezc zastoso- 
wanie narzedzia numeryczne opracowane na potrzeby niniejszej rozprawy 
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